%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Mach-Zehnder modulator" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: David Rörich
% date: 27.11.2013
%
% EXAMPLES: 
% (1) Transfer characteristic
% mach_zehnder_modulator([1 0], -0.5, 'push_pull', 1000, 1e5, 'single_tone', 100, 16, 1e9, 'rc', 'field', [0 0], 1, 0, 'dummy.png')
% (2) Spectrum of output field for single tone driving signal with all insets
% mach_zehnder_modulator([0 1], 0, 'push_pull', 1000, 1e5, 'single_tone', 100, 16, 1e9, 'rc', 'field', [1 1], 1, 0, 'dummy.png')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN main function
function mach_zehnder_modulator(...
    plot_en,...             % plot enable vector [transfer ch., spectrum, iqm]
    U_dc,...                % DC bias voltage, normalized to U_pi (halfwave switching voltage)
    operation_mode,...      % 'push_push' or 'push_pull'
    extinction_ratio_db,... % extinction ratio in dB
    num_samples,...         % number of samples in driving signals
    signal_type,...         % 'single_tone' or 'OOK'
    modulation_index,...    % in percent
    oversampling,...        % oversampling factor
    symbol_rate,...         % in GBd (or single tone freq in GHz for single_tone)
    pulse_shape,...         % 'rect' or 'rc'
    psd_mode,...            % 'field' or 'intensity'
    subplot_en,...          % inset enabling vector: [transfer_char time_domain]
    svg_en,...              % if 1, then also create svg file
    webdemo_en,...          % if 0, no file output but gnuplot display
    filename...             % name of output file
)

% Number of given paramters must exactly match the num of possible parameters
if nargin < nargin('mach_zehnder_modulator')
    error('Too few arguments.');
end
if nargin > nargin('mach_zehnder_modulator')
    error('Too many arguments.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters that are currently constant
U_pi = 3; % halfwave switching voltage in Volts
U_dc = U_dc * U_pi; % DC bias in Volts
input_power = 1e-3; % optical input power in Watt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Unit and dimension conversions
symbol_rate = symbol_rate * 1e9;

% Parameter checks
if (num_samples > 1e7)
    num_samples = 1e7;
    warning('"num_samples is limited to 1e7. Setting "num_samples" to 1e7 now.');
end
if (num_samples <= 0)
    num_samples = 1e4;
    warning('"num_samples" must be greater 0 and is now set to 1e4.');
end
if (symbol_rate <= 0)
    symbol_rate = 1e9;
    warning('"symbol_rate" must be greater 0 and is now set to 1e9.');
end
if (modulation_index <= 0)
    modulation_index = 100;
    warning('"modulation_index" must be greater 0 and is now set to 100.');
end
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, set up graphical output

% Clear figure objects (useful when executed at octave prompt)
clf;

% Create an invisible figure.
if (webdemo_en)
    fig = figure(1); set(fig, 'visible', 'off');
else
    fig = figure(1); set(fig, 'visible', 'on');
end

landscape = '-S930,350';
portrait  = '-S640,480';

output_format = landscape;     % default

set(0, 'defaultlinelinewidth', 2);

if (webdemo_en)
    set(0, 'defaultaxesfontsize', 8);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 8);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 2;
elseif (svg_en)
    set(0, 'defaultaxesfontsize', 6);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 6);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 1;
else
    set(0, 'defaultaxesfontsize', 10);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 8);
    set(0, 'defaulttextfontname', 'Helvetica');
    const_lw = 2;
end

% define colors
uni_mittelblau = [0 126/255 198/255];
uni_rot = [238/255 28/255 35/255];
uni_apfelgruen = [140/255 198/255 62/255];

% END, set up graphical output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plot transfer characteristics
if plot_en(1)
    plot_transfer_characteristics(num_samples, input_power, U_pi, U_dc,...
        operation_mode, extinction_ratio_db, [uni_mittelblau; uni_rot]);

    % Set aspect ratio of image file
    output_format = landscape;
end

% Plot spectrum of modulated signal
if plot_en(2)
    plot_mzm_output_spectrum(signal_type, symbol_rate, num_samples,...
    pulse_shape, oversampling, modulation_index, extinction_ratio_db,...
    input_power, [uni_mittelblau; uni_rot; 0.5*uni_apfelgruen], psd_mode,...
    subplot_en);
    
    % Set aspect ratio of image file
    output_format = portrait;
end

% Create image file that will be displayed on the webpage. Its size
% must be 640x480 for "normal" and 1280x960 for wide demo slides.
if(webdemo_en)
    print(filename, '-dpng',output_format);
end
if (svg_en)
    print(['fig_',filename], '-dsvg',output_format);
end

end
% END main function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN Function definitions
function E_out = mzm_modulate( u_1, u_2, E_0, U_dc_1, U_dc_2, extinction_ratio_db )
%MZM_MODULATE Modulate an optical carrier using a Mach-Zehnder modulator.
%
%INPUT:
%  u_1 Driving voltage of upper phase modulator, normalized to voltage U_pi
%      that causes phase shift of pi.
%  u_2 Driving voltage of lower phase modulator, normalized to voltage U_pi
%      that causes phase shift of pi.
%  E_0 Amplitude of electrical input field of optical carrier.
%  U_dc_1 Constant bias voltage of upper phase modulator, normalized to 
%      voltage U_pi that causes phase shift of pi.
%  U_dc_2 Constant bias voltage of lower phase modulator, normalized to 
%      voltage U_pi that causes phase shift of pi.
%  extinction_ratio_db Extinction Ratio in dB.
%
%OUTPUT:
%  E_out Complex envelope of output electrical field.

% Splitting factors
if (abs(extinction_ratio_db) <= 1000)
    er = 10^(extinction_ratio_db / 10);
    epsilon = 0.5*sqrt(1 - ((er - 1) / (er + 1))^2);
else % prevent overflow of er
    epsilon = 0.0;
end
Y_1 = sqrt(0.5 + epsilon); % upper arm splitting factor
Y_2 = sqrt(0.5 - epsilon); % lower arm splitting factor

% phase upper arm
phi_1 = (u_1 + U_dc_1) * pi;
% phase lower arm
phi_2 = (u_2 + U_dc_2) * pi;

% Complex envelope of output electrical field
E_out = E_0 / sqrt(2) * (Y_1*exp(1i*phi_1) + Y_2*exp(1i*phi_2));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plot_transfer_characteristics(...
    num_samples,...         % number of samples in driving signals
    input_power,...         % power of optical carrier
    U_pi,...                % halfwave switching voltage in Volts
    U_dc,...                % DC bias voltage in Volts
    operation_mode,...      % 'push_push' or 'push_pull'
    extinction_ratio_db,... % Extinction ratio in dB
    line_colors )           % Line colors for plots
% PLOT_TRANSFER_CHARACTERISTICS Plot field, power and phase transfer 
% characteristic of a Mach-Zehnder modulator.
%
% The driving voltage u_upper of the upper arm of the MZM is u_1 + U_dc with
% u_1 = -U_pi ... U_pi. The driving voltage u_lower of the lower arm of the 
% MZM depends on the operation mode:
%   In push-push mode: u_lower = u_upper
%   In push-pull mode: u_lower = -u_upper

% Input electrical field amplitude
E_0 = sqrt(2*input_power);

% driving voltage upper arm (in Volts)
u_1 = linspace(-1, 1, num_samples) * U_pi;

% DC bias in upper arm (in Volts)
U_dc_1 = U_dc;

% driving voltage lower arm (in Volts)    
if strcmp(operation_mode, 'push_pull')
    u_2 = -u_1;
    U_dc_2 = - U_dc_1;
elseif strcmp(operation_mode, 'push_push')
    u_2 = u_1;
    U_dc_2 = U_dc_1;
else
    error(['Unkwnown operation mode "' operation_mode...
        '". Possible operation modes: "push_pull" and "push_push"']);
end

% Complex envelope of output electrical field
E_out = mzm_modulate(u_1/U_pi, u_2/U_pi, E_0, U_dc_1/U_pi, U_dc_2/U_pi,...
    extinction_ratio_db);

% PLOTS
set(gcf, 'position', [0 0 930 350]);
% Plot the power
subplot(1, 2, 1);
plot(u_1/U_pi, 0.5 * abs(E_out).^2 / input_power, 'color', line_colors(1,:));
hold all;    
% Plot the amplitude. The amplitude is defined to be the negative magnitude
% where the cosine of the angle is negative.
amplitude_sign = ones(size(E_out));
amplitude_sign( cos(angle(E_out)) >= 0 ) = 1;
amplitude_sign( cos(angle(E_out)) < 0 ) = -1;
plot(u_1/U_pi, amplitude_sign .* abs(E_out) / E_0, 'color', line_colors(2,:));
xlabel('u_1/U_{\pi}');
title('Power and field transfer characteristic');
legend('P_{out}/P_{in}', 'E_{out}/E_{in}', 'location', 'southeast');
hold off;
grid on;
ylim([-1.1, 1.1]);    
set(gca, 'position', [0.05 0.15 0.4 0.72]);
% Plot the phase
phase_trans = unwrap(angle(E_out))/pi;
phase_range = max(phase_trans) - min(phase_trans);
subplot(1, 2, 2);
plot(u_1/U_pi, phase_trans, 'color', line_colors(1,:));
title('Phase characteristic');
xlabel('u_1/U_{\pi}');
ylabel('Phase/\pi in rad');
set(gca, 'YLim', [min(phase_trans)-0.05*phase_range max(phase_trans)+0.05*phase_range]);
grid on;    
set(gca, 'position', [0.57 0.15 0.4 0.72]);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function excess_loss = plot_mzm_output_spectrum(...
    signal_type,...         % 'OOK' or 'single_tone'
    symbol_rate,...         % OOK: symbols/s, single tone: oscillation freq. in Hz
    num_samples,...         % number of samples in the driving signal
    pulse_shape,...         % OOK: 'rect' or 'rc' (raised cosine), single tone: ignored
    oversampling,...        % OOK: pulse width in samples, single tone: ignored
    modulation_index,...    % Ratio of max signal voltage to full driving voltage in percent
    extinction_ratio_db,... % MZM extinction ratio in dB
    optical_input_power,... % Power of input optical carrier
    line_colors,...         % Vector of line colors for plot
    psd_mode,...            % 'field' or 'intensity'
    subplot_en )            % inset enabling vector: [transfer_char time_domain]
% PLOT_MZM_OUTPUT_SPECTRUM Plot spectrum of modulated optical carrier and
% calculate the modulation excess loss in dB.
% 
% Plots the power spectral density (PSD) of the driving signal and optical
% output of a Mach-Zehnder modulator driven by a OOK or single tone signal. If
% the signal type is OOK, the modulator is driven in the minimum transmission
% point. If the signal type is single tone it is driven in the quadrature
% point. 
%
% If psd_mode is 'field' the PSD of the output electrical field is
% plottet. If psd_mode is 'intensity' the PSD of the output light intensity
% is plotted (square of electrical field).
%
% The modulation index is defined as the ratio of the maximum driving voltage 
% U_max to the voltage U_f that drives the MZM to full scale. U_f and U_max 
% are defined as follows:
% OOK - minimum transmission point - U_f = U_pi/4
%                                    U_max = 0.25 * U_pi * modulation_index / 100
% single tone - quadrature point   - U_f = U_pi/2
%                                    U_max = 0.5 * U_pi * modulation_index / 100
% For both signal types the voltage swing is -U_max ... U_max.

% Amplitude of optical input field
E_0 = sqrt(2*optical_input_power);

% Inset enable flags
trans_inset_en = subplot_en(1); % transfer char.
td_inset_en = subplot_en(2); % time domain signals

% Positions of insets
if trans_inset_en && td_inset_en;
    trans_inset_pos = [0.47 0.5 0.2 0.2];
    trans_rect_pos = [0.39 0.45 0.27 0.3];
else
    trans_inset_pos = [0.72 0.5 0.2 0.2];
    trans_rect_pos = [0.65 0.45 0.24 0.3];
end
drive_inset_pos = [0.47 0.23 0.2 0.2];
drive_rect_pos = [0.39 0.175 0.27 0.3];
out_inset_pos = [0.72 0.5 0.2 0.2];
out_rect_pos = [0.66 0.45 0.23 0.3];
rect_face_color = 0.95*[1 1 1];

% Line color of driving signal
drive_line_color = line_colors(3,:);

% Some variables that depend on psd_mode
if (strcmp(psd_mode, 'field'))
    plot_field_psd = 1;
    out_legend_str = 'output electrical field';
    out_line_color = line_colors(2,:);
    
elseif (strcmp(psd_mode, 'intensity'))
    plot_field_psd = 0;
    out_legend_str = 'output intensity';
    out_line_color = line_colors(1,:);
else
    error(['Parameter "psd_mode" must be "field" or "intensity", found ',...
        psd_mode]);
end

if(strcmp(signal_type, 'single_tone')) % Modulate with single tone
    % Amplitude of driving signal
    A = 0.5 * (modulation_index / 100);
    
    single_tone_freq = symbol_rate;
    sampling_freq = oversampling * single_tone_freq; % in GHz
    t = linspace(0, num_samples/sampling_freq, num_samples);
    
    u_1 = A * cos(2*pi*single_tone_freq*t);

    % Operate in quadrature point
    E_out = mzm_modulate(u_1, -u_1, E_0, -0.5, 0.5, extinction_ratio_db);
    
    % Generate transfer characteristic for reference.
    u_trans = A * linspace(-1, 1, num_samples);    
    trans = mzm_modulate(u_trans, -u_trans, E_0, -0.5, 0.5, extinction_ratio_db);
    if plot_field_psd
        % field transfer charactersitic
        amplitude_sign = ones(size(trans));
        amplitude_sign( cos(angle(trans)) >= 0 ) = 1;
        amplitude_sign( cos(angle(trans)) < 0 ) = -1;
        trans = amplitude_sign .* abs(trans) / E_0;
        subplot_title = 'Field transfer char.';
    else
        % power transfer characteristic
        trans = 0.5 * abs(trans).^2 / optical_input_power;
        subplot_title = 'Power transfer char.';
    end
    
    % Calculate estimated PSDs
    [P_uu, w_u] = spectrum(u_1);
    if plot_field_psd
        [P_EE, w_E] = spectrum(E_out);
    else
        [P_EE, w_E] = spectrum(E_out .* conj(E_out));
    end

    % Find power of carrier in each PSD    
    P_carrier_u = max(P_uu( w_u*sampling_freq/2/pi > 0.9 * single_tone_freq &...
        w_u*sampling_freq/2/pi < 1.1 * single_tone_freq));
    P_carrier_E = 1;
    if plot_field_psd
        P_carrier_E = max(P_EE( w_E*sampling_freq/2/pi > 0.9 * single_tone_freq &...
            w_E*sampling_freq/2/pi < 1.1 * single_tone_freq));
    else
        % When plotting intensity, carrier freq. is doubled due to squaring
        P_carrier_E = max(P_EE( w_E*sampling_freq/2/pi > 0.9 * 2 *single_tone_freq &...
            w_E*sampling_freq/2/pi < 1.1 * 2 * single_tone_freq));
    end
        
    % Calculate mean power of optical output signal
    P_E_mean = sum(abs(E_out).^2) / length(E_out);
    
    % Calculate modulation excess loss
    excess_loss = 10*log10(optical_input_power/(P_E_mean/2));

    % Plot    
    plot(w_E*sampling_freq/2/pi, 10*log10(P_EE/P_carrier_E), 'color', ...
        out_line_color);
    hold all;
    plot(w_u*sampling_freq/2/pi, 10*log10(P_uu/P_carrier_u), 'color', ...
        drive_line_color, 'linestyle', '--', 'marker', 'none', 'markerfacecolor', 'none');
    hold off;
    title(['Power spectral density. Excess loss: ',...
        sprintf('%.1f dB', excess_loss)]);
    ylabel('Power spectral density in dBc');
    xlabel('Frequency in Hz');
    ylim([-100 0]);
    legend(out_legend_str, 'driving signal');
    grid('on');
    grid('minor', 'off');
    main_ax = gca;    
    if trans_inset_en        
        sub_ax = axes('position', trans_inset_pos);
        plot(u_trans, trans, 'color', out_line_color);
        set(sub_ax, 'color', 'w');        
        set(sub_ax, 'XTick', [-A 0 A]);
        set(sub_ax, 'XLim', [-1.1*A 1.1*A]);
        set(sub_ax, 'XTickLabel', {'u_{min}'; '0'; 'u_{max}'});
        set(sub_ax, 'YLim', [-1.1 1.1]);
        set(sub_ax, 'YTick', [-1 0 1]);
        set(sub_ax, 'XGrid', 'on');
        set(sub_ax, 'YGrid', 'on');
        set(sub_ax, 'Title', subplot_title);      
    end
    
elseif(strcmp(signal_type, 'OOK'))

    sampling_freq = symbol_rate * oversampling;

    %% OOK signal and spectrum
    num_symbols = ceil(num_samples / oversampling);
    bits = randi(2, 1, num_symbols) - 1;
    % Upsampling
    bits_up = zeros(1, length(bits)*oversampling);
    bits_up(1:oversampling:end) = bits;
    if(strcmp(pulse_shape, 'rect'))
        h_pulse = ones(1, oversampling); % rectangular
    elseif(strcmp(pulse_shape, 'rc'))
        % raised cosine impulse with roll-off factor 0.5 and 50% overlap
        h_pulse = rc_impulse_response(oversampling*2+1,0.5);
    else
        error(['Unknown pulse shape. Supported: "rect", "rc". Found: '...
            pulse_shape]);
    end
    pulses = conv(bits_up, h_pulse);
    pulses = 0.5 * pulses - 0.25;
    u_1 = pulses * modulation_index / 100;

    % Operate in minimum transmission point
    E_out = mzm_modulate(u_1, -u_1, E_0, -0.25, 0.25, extinction_ratio_db);    

    % Generate transfer characteristic for reference.
    u_trans = linspace(min(u_1), max(u_1), num_samples);    
    trans = mzm_modulate(u_trans, -u_trans, E_0, -0.25, 0.25, extinction_ratio_db);
    if plot_field_psd
        % field transfer characteristic
        amplitude_sign = ones(size(trans));
        amplitude_sign( cos(angle(trans)) >= 0 ) = 1;
        amplitude_sign( cos(angle(trans)) < 0 ) = -1;
        trans = amplitude_sign .* abs(trans) / E_0;
        subplot_title = 'Field transfer char.';
    else
        % power transfer characteristic
        trans = 0.5 * abs(trans).^2 / optical_input_power;
        subplot_title = 'Power transfer char.';
    end
    
    % Calculate estimated PSDs
    [P_uu, w_u] = spectrum(u_1);
    if plot_field_psd
        [P_EE, w_E] = spectrum(E_out);
    else
        [P_EE, w_E] = spectrum(E_out .* conj(E_out));
    end

    % Calculate modulation excess loss
    P_E_mean = sum(abs(E_out).^2) / length(E_out);
    excess_loss = 10*log10(optical_input_power/(P_E_mean/2));
    
    % Normalization of PSDs (removing physical units from PSDs)
    PSD_u_norm = mean(P_uu);
    PSD_E_norm = mean(P_EE);
    
    % Plot
    plot(w_E*sampling_freq/2/pi, 10*log10(P_EE/PSD_E_norm), 'color',...
        out_line_color);
    hold on;
    plot(w_u*sampling_freq/2/pi, 10*log10(P_uu/PSD_u_norm), 'color',...
        drive_line_color);
    hold off;
    title(['Power spectral density (w.r.t. mean PSD). Excess loss: ',...
        sprintf('%.1f dB', excess_loss)]);
    legend(out_legend_str, 'driving signal');
    ylabel('Power spectral density in dB');
    xlabel('Frequency in Hz');
    ylim([-50 50]);
    xlim([-symbol_rate/10 max(w_u*sampling_freq/2/pi)]); % create a small space at left margin of plot
    grid('on');
    grid('minor', 'off');
    main_ax = gca;
    if trans_inset_en
        sub_ax = axes('position', trans_inset_pos);
        plot(u_trans, trans, 'color', out_line_color);
        set(sub_ax, 'color', 'w');
        set(sub_ax, 'XTick', [min(u_1) 0 max(u_1)]);
        set(sub_ax, 'XLim', [1.1*min(u_1) 1.1*max(u_1)]);
        set(sub_ax, 'XTickLabel', {'u_{min}'; '0'; 'u_{max}'});
        set(sub_ax, 'YLim', [-0.1 1.1]);
        set(sub_ax, 'YTick', [0 0.5 1]);
        set(sub_ax, 'XGrid', 'on');
        set(sub_ax, 'YGrid', 'on');
        set(sub_ax, 'Title', subplot_title);
    end
else
    error(['Unknown signal type. Supported: "single_tone" and "OOK". Found: ' ...
        signal_type]);
end

% Plot time domain signals in two insets
if td_inset_en    
    out_ax = axes('position', out_inset_pos);
    M = 5*oversampling; % no of samples to plot
    t_vec = (0:M-1)/(oversampling*sampling_freq);
    if plot_field_psd
        % plot output field        
        amplitude_sign = ones(size(E_out(1:M)));
        amplitude_sign( cos(angle(E_out(1:M))) >= 0 ) = 1;
        amplitude_sign( cos(angle(E_out(1:M))) < 0 ) = -1;
        out_sig = amplitude_sign .* abs(E_out(1:M))/E_0;
        plot(t_vec, out_sig, 'color', out_line_color);
    else
        % plot output power
        out_sig = 0.5 * abs(E_out(1:M)).^2 / optical_input_power;
        plot(t_vec, out_sig, 'color', out_line_color);
    end
    if (strcmp(signal_type,'single_tone'))
        out_ylim = [-1.1 1.1];
        out_y_ticks = [-1 0 1];
    else
        out_ylim = [-0.1 1.1];
        out_y_ticks = [0 0.5 1]
    end
    set(out_ax, 'color', 'w');
    set(out_ax, 'XTick', [t_vec(1) t_vec(oversampling) t_vec(4*oversampling)]);
    set(out_ax, 'XLim', [t_vec(1) t_vec(end)]);
    set(out_ax, 'XTickLabel', [{'0', 'T_S', '4T_S'}]);
    set(out_ax, 'YLim', out_ylim);
    set(out_ax, 'YTick', out_y_ticks);
    set(out_ax, 'YTickLabel', []);
    set(out_ax, 'XGrid', 'on');
    set(out_ax, 'YGrid', 'on');
    set(out_ax, 'Title', out_legend_str);
    
    in_ax = axes('position', drive_inset_pos);
    plot(u_1(1:M), t_vec, 'color', drive_line_color);
    set(in_ax, 'color', 'w');
    set(in_ax, 'XTick', [min(u_1) 0 max(u_1)]);
    set(in_ax, 'XLim', [1.1*min(u_1) 1.1*max(u_1)]);
    set(in_ax, 'XTickLabel', {'u_{min}'; '0'; 'u_{max}'});
    set(in_ax, 'YLim', [t_vec(1) t_vec(end)]);
    set(in_ax, 'YTick', [t_vec(1) t_vec(oversampling) t_vec(4*oversampling)]);
    set(in_ax, 'YTickLabel', [{'0', 'T_S', '4T_S'}]);
    set(in_ax, 'XGrid', 'on');
    set(in_ax, 'YGrid', 'on');
    set(in_ax, 'Title', '');
    set(in_ax, 'YLabel', 'time');
    set(in_ax, 'YAxisLocation', 'left');
end
% Plot rectangles as background of inset axes.
axes(main_ax);
if trans_inset_en
    rectangle('Position', fig2ds_pos(main_ax, trans_rect_pos), 'FaceColor', ...
        rect_face_color, 'EdgeColor', 'none');
end
if td_inset_en
    rectangle('Position', fig2ds_pos(main_ax, drive_rect_pos), 'FaceColor', ...
        rect_face_color, 'EdgeColor', 'none');
    rectangle('Position', fig2ds_pos(main_ax, out_rect_pos), 'FaceColor', ...
        rect_face_color, 'EdgeColor', 'none');
end

end % function end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function h = rc_impulse_response(N, a)
% RC_IMPULSE_RESPONSE Generate a raised-cosine shaped impulse.
%
% h = RC_IMPULSE_RESPONSE(N, a) returns a raised-cosine impulse h with roll-off 
% factor a and length N. Note that this is not the impulse response of a raised-
% cosine filter but an impulse that has the same shape as the transfer function 
% of a raised-cosine filter.
h = zeros(1, N); % init
t_n_vec = linspace(-1, 1, N); % discrete time
for k = 1:N
    t_n = t_n_vec(k);
    if (abs(t_n) <= (1-a)/2)
        h(k) = 1;
    elseif( ((1 - a)/2) < abs(t_n) && abs(t_n) <= (1 + a)/2)
        h(k) = 0.5 * (1 + cos(pi/a * (abs(t_n) - (1 - a)/2)));
    else
        h(k) = 0;
    end
end
end

function [P_xx, w] = spectrum(x)
% SPECTRUM
%
% [P_xx, w] = SPECTRUM(x) returns the one-sided spectral power density of x, 
% P_xx, and angular frequency vector w. 
    % x must be column vector
    if size(x,2) > size(x,1)
        x = x';
    end
    
    % Calculate one-sided spectrum. x is weighted with Hamming window.
    [P_xx, w] = periodogram(x, hamming(length(x)), 'onesided');
    
    % Smoothen spectrum by taking the average of every chunk_size freq bins.
    % If P_xx is not an integer multiple of chunk_size, the remaining bins
    % are discarded.
    chunk_size = 100;
    num_chunks = floor(length(P_xx)/chunk_size);
    P_xx_smooth = zeros(num_chunks, 1);
    w_smooth = zeros(num_chunks, 1);
    for k = 1:num_chunks
        P_xx_smooth(k) = mean( P_xx( (k-1)*chunk_size+1:k*chunk_size ) );
        w_smooth(k) = w((k-1)*chunk_size + ceil(chunk_size/2));
    end
    P_xx = P_xx_smooth;
    w = w_smooth;
end

function ds_pos = fig2ds_pos(ah, pos)
% FIG2DS_POS Convert a position from normalized figure to data space units.
%
% IN:
%   ah  Handle of axes that contains the dataspace
%   pos Four element position vector [left bottom width height] in normalized
%       figure units. (Values between 0 and 1.0).
% OUT:
%  ds_pos Four element position vector in data space units of the given axes.

    % Position that is to be converted
    x = pos(1);
    y = pos(2);
    width = pos(3);
    height = pos(4);
    
    % Limits of data space of given axes
    ax_xlim = get(ah, 'XLim');
    ax_ylim = get(ah, 'YLim');
    dsx_max = ax_xlim(2);
    dsx_min = ax_xlim(1);
    dsy_max = ax_ylim(2);
    dsy_min = ax_ylim(1);
    
    % Position of given axes inside the figure
    ax_pos = get(ah, 'Position');
    fig_x = ax_pos(1);
    fig_y = ax_pos(2);
    fig_width = ax_pos(3);
    fig_height = ax_pos(4);
    
    % Conversion
    ds_pos(1) = (dsx_max - dsx_min)/fig_width * (x - fig_x) + dsx_min;
    ds_pos(2) = (dsy_max - dsy_min)/fig_height * (y - fig_y) + dsy_min;
    ds_pos(3) = (dsx_max - dsx_min)/fig_width * width;
    ds_pos(4) = (dsy_max - dsy_min)/fig_height * height;
end

% Not used.
function plot_axes_margin(ax_parent, ax_child_pos, margin, face_color)
    if nargin < 4
        face_color = [1 1 1];
    end
    % position of rectangle in normalized axes units
    x_ax = ax_child_pos(1) - margin;
    y_ax = ax_child_pos(2) - margin;
    width_ax = ax_child_pos(3) + 2*margin;
    height_ax = ax_child_pos(4) + 2*margin;
    rect_pos_ax = [x_ax y_ax width_ax height_ax];
    disp(rect_pos_ax);
    
    % convert to dataspace units
    rect_pos_ds = fig2ds_pos(ax_parent, rect_pos_ax);
    disp(rect_pos_ds);
    
    % ensure that the margin rectangle does not exceed the limits of the parent axes
    ax_xlim = get(ax_parent, 'XLim');
    ax_ylim = get(ax_parent, 'YLim');
    disp(ax_xlim);
    disp(ax_ylim);
    x_max = ax_xlim(2);
    x_min = ax_xlim(1);
    y_max = ax_ylim(2);
    y_min = ax_ylim(1);
    rect_pos_ds(1) = min(rect_pos_ds(1), x_max); % x1 not larger than x_max
    rect_pos_ds(1) = max(rect_pos_ds(1), x_min); % x1 not less than x_min
    rect_pos_ds(2) = min(rect_pos_ds(2), y_max); % y1 not larger than y_max
    rect_pos_ds(2) = max(rect_pos_ds(2), y_min); % y1 not less than y_min
    rect_pos_ds(3) = min(rect_pos_ds(3), x_max - rect_pos_ds(1)); % x2 not larger than x_max
    rect_pos_ds(3) = max(rect_pos_ds(3), x_min - rect_pos_ds(1)); % x2 not less than x_min
    rect_pos_ds(4) = min(rect_pos_ds(4), y_max - rect_pos_ds(2)); % y2 not larger than y_max
    rect_pos_ds(4) = max(rect_pos_ds(4), y_min - rect_pos_ds(2)); % y2 not less than y_min
    
    disp(rect_pos_ds);
    % plot rectangle
    rectangle('Position', rect_pos_ds, 'FaceColor', face_color,...
            'EdgeColor', 'none');
end
% END Function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
