%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN main program (control)
function [min, avg] = spectral_shaping(...
    M,              ...% size of QAM/PAM alphabet
    qamMode,        ...% mode of QAM/PAM: 'bpsk', 'qpsk', 'qam', 'pam'
    Nsub,           ...% number of subcarriers for OFDM
    Nsym,           ...% number of OFDM-symbols to be generated
    fcz,            ...% first subcarrier to cancel
    lcz,            ...% last subcarrier to cancel
    ecz,            ...% number of edge subcarriers to disable
    ovsamp,         ...% factor of oversampling for aic-calculation
    Nup,            ...% factor of oversampling used for finer plot of spectrum
    mode,           ...% mode of the cancellation:
                    ...% 'none', 'dt'/'tn', 'aiconly', 'aic', 'fir', 'filtercmp'
    Nc,             ...% number of subcarriers used for cancellation signal
    fir_params,     ...% parameters for FIR filter design:
                    ...% [passband_db stopband_db nb_taps transition_width neg_en]
    iq_imb,         ...% IQ imbalance, set to 0 to disable, otherwise pass a 
                    ...% vector [phi_imb a_imb dc_imb], a_imb and dc_imb are optional
    plot_en,        ...% enable plotting [primary secondary cmp_pm cmp_ls cmp_rc]
    dB_en,          ...% if 1, then plot in dB
    webdemo_en,     ...% if 0, then test-mode in Octave/Matlab commandline
    filename,       ...% output file for webdemo
    filetype,       ...% Type of output image: 'png', 'eps', 'svg'
    yscale,         ...% optional ymin and ymax
    xscale,         ...% optional xmin and xmax
    heading         ...% optional title
  )
  
  % set number of required parameters
  nargin_default = 18;
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN check input parameters for consitency  
  
  % start measurement for timings of algorithm parts
  tic
  
  % convert negative subcarriers in range from -Nsub/2 to -1 
  % to positive subcarriers in range from Nsub/2 to Nsub-1
  fcz = mod(fcz, Nsub);
  lcz = mod(lcz, Nsub);
  
  % sanity check of carrier deactivation
  if fcz > Nsub || lcz > Nsub || 2*ecz >= Nsub
    fprintf('Error: please check selected/deactivated carriers\n');
    return
  end
    
  % set default fir filter implementation
  if strcmpi(mode,'fir')
    mode = 'firls';
  end

  % END check input parameters for consitency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN other initializations
  
  % set true to enable logging
  log_en = ~webdemo_en;
  
  % compatibility mode to run this script in Matlab
  matlab_mode_en = false;
  
  % set language for plot labels, choose between 'de' and 'en'
  language = 'en';
    
  % if enabled little dots will be drawn on the bottom of the plot 
  % to indicate disabled subcarriers
  carrier_info_plot_en = true;
  
  % set true to append nb of taps to filename
  fir_length_filename_en = false;

  % calculates the minimum and mean notch-depth of the OFDM signal
  power_calc_en = false;
  
  % reset random generator for caching if webdemo enabled
  if webdemo_en
    if matlab_mode_en
      % alternative method when using Matlab instead of Octave
      rng(0353098645);
    else 
      rand('state',1:500);
    end
  end
  
  % convert QPSK and BPSK inputs to QAM and PAM forms
  switch lower(qamMode)
    case 'qpsk'
      M = 4;
      qamMode = 'qam';
    case 'bpsk'
      M = 2;
      qamMode = 'pam';
  end
  
  % expand plot vectors if necessary
  if length(plot_en) == 1
    plot_en = [plot_en 0 0 0 0];
  elseif length(plot_en) == 2
    plot_en = [plot_en 0 0 0];
  end
  
  % prepare iq imbalance params
  switch length(iq_imb)
    case 1
      if iq_imb == 0
        imb_en  = 0;
      else
        imb_en  = 1;
        phi_imb = iq_imb(1);
        a_imb   = 1;
        dc_imb  = 0;
      end
    case 2
      imb_en  = 1;
      phi_imb = iq_imb(1);
      a_imb   = iq_imb(2);
      dc_imb  = 0;
    case 3
      imb_en  = 1;
      phi_imb = iq_imb(1);
      a_imb   = iq_imb(2);
      dc_imb  = iq_imb(3);
  end
  
  % default values, may be overwritten for certain cancellation modes
  opt.disableZeroCarrier = true;
  opt.disableEdgeCarriers = ecz ~= 0;
  
  opt.disableNotchCarriers = false;
  opt.disableNotchCarriersNeg = false;
  opt.activateAIC = false;
  opt.onlyAIC = false;
  opt.applyFIR = false;
  opt.filtercmp_en = false;
  
  % switching of cancellation mode
  switch lower(mode)
    case {'firpm', 'firls', 'firrc'}
      opt.applyFIR = true;
      opt.disableNotchCarriers = true;
      if length(fir_params) >= 5 && fir_params(5) == 1
        opt.disableNotchCarriersNeg = true;
      end
    case 'aic'
      opt.activateAIC = true;
      opt.disableNotchCarriers = true;
    case {'dt','disabledtones','tn','tonenulling'}
      opt.disableNotchCarriers = true;
    case {'aiconly', 'onlyaic'}
      opt.activateAIC = true;
      opt.onlyAIC = true;
    case {'filtercmp'}
      plot_en(1) = 0;
      plot_en(2) = 0;
      opt.filtercmp_en = true;
      carrier_info_plot_en = false;
  end
  
  % END other initializations
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN set up graphical output
  
  if sum(plot_en)
    if webdemo_en
      % create an invisible figure
      fig = figure(1); set(fig, 'visible', 'off');
    else
      fig = figure(1); set(fig, 'visible', 'on');
    end
    
    output_format  = '-S640,480';
    
    legendFontSize = 9;
    
    set(fig, 'defaultlinelinewidth', 2);
    switch lower(filetype)
      case 'png'
        if webdemo_en
          set(fig, 'defaultaxesfontsize', 8);
          set(fig, 'defaultaxesfontname', 'Helvetica');
          set(fig, 'defaulttextfontsize', 8);
          set(fig, 'defaulttextfontname', 'Helvetica');
        else
          set(fig, 'defaultaxesfontsize', 10);
          set(fig, 'defaultaxesfontname', 'Helvetica');
          set(fig, 'defaulttextfontsize', 8);
          set(fig, 'defaulttextfontname', 'Helvetica');
        end
      case 'svg'
        W = 14;
        H = W*0.75;
        set(fig,'PaperUnits','centimeters');
        set(fig,'PaperOrientation','portrait');
        set(fig,'PaperSize',[H,W]);
        set(fig,'PaperPosition',[0,0,W,H]);
        
        set(fig, 'defaultaxesfontsize', 6);
        set(fig, 'defaultaxesfontname', 'Times');
        set(fig, 'defaulttextfontsize', 6);
        set(fig, 'defaulttextfontname', 'Times');
        set(fig, 'defaultlinelinewidth', 1);
      case {'eps', 'pdf'}
        W = 14;
        H = W*0.75;
        set(fig,'PaperUnits','centimeters');
        set(fig,'PaperOrientation','portrait');
        set(fig,'PaperSize',[H,W]);
        set(fig,'PaperPosition',[0,0,W,H]);
        
        fontSize = 13;
        legendFontSize = 13;
        set(fig, 'defaultaxesfontsize', fontSize);
        set(fig, 'defaultaxesfontname', 'Helvetica');
        set(fig, 'defaulttextfontsize', fontSize);
        set(fig, 'defaulttextfontname', 'Helvetica');
        set(fig, 'defaultlinelinewidth', 1);
        
        if strcmpi(filetype, 'eps')
            filetype = 'epsc';
        end
        
        if matlab_mode_en
          box on
          ax = gca;
          ax.BoxStyle = 'full';
        end
    end
    
    % define colors for plot
    COLOR.BLUE  = [  0 126 198]/255;
    COLOR.RED   = [238  28  35]/255;
    COLOR.GREEN = [140 198  62]/255;
    COLOR.PINK  = [237   0 140]/255;
    COLOR.GREY  = [224 224 224]/255;
  end
  
  % END set up graphical output
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN calculation
  
  if plot_en(1)
    % generate Nsym vectors of Nsub random QAM-symbols
    console(log_en, 'generate QAM-symbols');
    s = generate_symbolvectors(Nsym, Nsub, M, qamMode);
    % calculation of cancellation matrix 
    if opt.activateAIC 
      console(log_en, 'caluclate AIC-matirx');
      C = calculate_aic_matrix(fcz, lcz, Nsub, Nc, ovsamp);
    else
      C = [];
    end
    
    % calculation of OFDM-signal
    console(log_en, 'calculationg OFDM-signal');
    ofdm_signal = calculate_OFDM_signal(s, Nsym, ecz, fcz, lcz, Nsub, Nc, opt, C);
    
    % generate and apply FIR-filter
    if opt.applyFIR
      console(log_en,'calculating FIR-filter');
      h = generate_fir_filter(Nsub, fcz, lcz, fir_params, mode);
      filter_len = length(h);
      if fir_length_filename_en 
        filename = [filename sprintf('_N=%i',filter_len)];
      end
      console(log_en,['apply FIR-filter of length N=' int2str(filter_len)]);
      ofdm_signal = conv(ofdm_signal, h);
    end
    
    % apply IQ-imbalance
    if imb_en
      ofdm_signal = apply_iq_imb(ofdm_signal, phi_imb, a_imb, dc_imb);
    end
    
    % calculation of spectrum
    console(log_en,'calculating OFDM spectrum');
    [ofdm_spectrum, ofdm_f_axis] = calculate_spectrum(ofdm_signal, Nsub, Nup, dB_en);
    
    % calculation for additional plots
    if plot_en(2)
      console(log_en,'calculation additional spectrum');
      if opt.activateAIC
        opt.onlyAIC = true;
        aic_signal = calculate_OFDM_signal(s, Nsym, ecz, fcz, lcz, Nsub, Nc, opt, C);
        [aic_spectrum, aic_f_axis] = calculate_spectrum(aic_signal, Nsub, Nup, dB_en);
      elseif opt.applyFIR
        [H, H_f_axis] = calculate_spectrum(h, Nsub, Nup, dB_en, 1);
      end
    end
    
  elseif opt.filtercmp_en
    if plot_en(3)
      pm_signal = generate_fir_filter(Nsub, fcz, lcz, fir_params, 'firpm');
      [pm_spectrum, pm_f_axis] = calculate_spectrum(pm_signal, Nsub, Nup, dB_en, 1);
    end
    if plot_en(4)
      ls_signal = generate_fir_filter(Nsub, fcz, lcz, fir_params, 'firls');
      [ls_spectrum, ls_f_axis] = calculate_spectrum(ls_signal, Nsub, Nup, dB_en, 1);
    end
    if plot_en(5)
      rc_signal = generate_fir_filter(Nsub, fcz, lcz, fir_params, 'firrc');
      [rc_spectrum, rc_f_axis] = calculate_spectrum(rc_signal, Nsub, Nup, dB_en, 1);
    end
  end
  
  % END calculation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN generate plots
  
  if sum(plot_en)
    console(log_en,'plotting');
    clf
    hold on
    grid on
    if plot_en(1)
      plot(ofdm_f_axis, ofdm_spectrum, 'Color', COLOR.BLUE);
      if plot_en(2)
        if opt.activateAIC
          plot(aic_f_axis, aic_spectrum, 'Color', COLOR.RED);
        elseif opt.applyFIR
          plot(H_f_axis, H, 'Color', COLOR.RED);
        end
      end
      lh = legend(get_legend(language, opt, plot_en(2)));
    elseif opt.filtercmp_en
      legend_strings = {};
      if plot_en(3)
        plot(pm_f_axis, pm_spectrum, 'Color', COLOR.BLUE);
        legend_strings = [legend_strings; {get_name(language, 'firpm')}];
      end
      if plot_en(4)
        plot(ls_f_axis, ls_spectrum, 'Color', COLOR.RED);
        legend_strings = [legend_strings;  {get_name(language, 'firls')}];
      end
      if plot_en(5)
        plot(rc_f_axis, rc_spectrum, 'Color', COLOR.GREEN);
        legend_strings = [legend_strings; {get_name(language, 'firrc')}];
      end
      lh = legend(legend_strings);
    end 
    
    set(lh, 'FontSize', legendFontSize);
    if plot_en(1)
        set(lh,'Location','southeast');
    end
    
    if strcmpi(mode, 'filtercmp')
        set(lh,'Location','northeast');
    end
        
    if power_calc_en        
      % calculate peak power:
      if dB_en
        tmpstr = 'dB';
      else
        tmpstr = '';
      end
      peak = max(ofdm_spectrum);
      %pre = blanks(3-ceil(log10((abs(peak))+1)));
      pre = '';
      console(log_en, [pre num2str(peak, '%.1f') tmpstr ' peak power']);    
      
      % calculate mean power in stopband:
      %f_stop = (Nup*fcz:Nup*lcz-1)-0.5*Nup*Nsub+1;
      f_stop = mod((Nup*fcz:Nup*lcz-1)-0.5*Nup*Nsub, Nup*Nsub)+1;
      X_sq_stop = ofdm_spectrum(f_stop);
      if dB_en
        X_sq_stop = 10.^(X_sq_stop./10);
      end
      min = max(X_sq_stop);
      avg = sum(X_sq_stop)/length(X_sq_stop);
      if dB_en
        min = 10*log10(min);
        avg = 10*log10(avg);
      end
      console(log_en, [num2str(min, '%.1f') tmpstr ' min notch depth']);
      console(log_en, [num2str(avg, '%.1f') tmpstr ' avg notch depth']);
    end
    
    % y axis limit and label 
    if nargin > nargin_default
      min_y = yscale(1);
      max_y = yscale(2);
    elseif plot_en(1)
      min_y = min(ofdm_spectrum);
      max_y = max(ofdm_spectrum);
      if dB_en
        min_y = min_y-5;
        max_y = max_y+5;
      else
        min_y = min_y-0.1;
        max_y = max_y+0.1;
      end
    else
      min_y = -60;
      max_y = +10;
    end
    
    if opt.filtercmp_en
      if dB_en
        ylabel('|H(f)|^2 in dB');
      else
        ylabel('|H(f)|');
      end
    else    
      if dB_en
        ylabel('|X(f)|^2 in dB');
      else
        ylabel('|X(f)|');
      end
    end
    
    ylim([min_y max_y]);
    
    % additional points for carrier information
    if carrier_info_plot_en
      if opt.disableEdgeCarriers
        draw_points(-Nsub/2:-Nsub/2+ecz, min_y, COLOR.BLUE);
        draw_points(+Nsub/2-ecz:+Nsub/2, min_y, COLOR.BLUE);
      end
      if opt.disableZeroCarrier
        draw_points(0, min_y, COLOR.BLUE);
      end
      if opt.activateAIC
         draw_points(to_Nsub_half(fcz-Nc:lcz+Nc, Nsub), min_y, COLOR.RED);
      elseif opt.disableNotchCarriers
        draw_points(to_Nsub_half(fcz:lcz,Nsub), min_y, COLOR.BLUE);
      end
      if opt.disableNotchCarriersNeg
        draw_points(to_Nsub_half(-lcz:-fcz, Nsub), min_y, COLOR.BLUE);
      end
    end
    
    % x axis limit and label
    min_x = -Nsub/2;
    max_x = +Nsub/2;
    
    if nargin > nargin_default+1
        min_x = xscale(1);
        max_x = xscale(2);
    end
    
    if nargin > nargin_default+2
      title(heading);
    end

    xlabel(get_name(language, 'xlabel'));
    xlim([min_x max_x]);
    if opt.filtercmp_en
      xtick = 1;
    else
      if rem(log2(Nsub), 1) == 0
        xtick = 2^(round(log2(Nsub/16)));
      elseif rem(Nsub, 100) == 0
        if Nsub < 801
          xtick = 50;
        else 
          xtick = 100;
        end
      else 
        xtick = 2^(round(log2(Nsub/16))); 
      end
    end
    set(gca,'XTick',floor(min_x/xtick)*xtick:xtick:ceil(max_x/xtick)*xtick);

    hold off
  end
  
  % END generate plots
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN output to file
  if sum(plot_en) && webdemo_en && ~strcmpi(filetype, 'none') 
      if strcmpi(filetype, 'png')
        print(filename, ['-d', filetype], output_format);
      else
        print(filename, ['-d', filetype], '-loose');
      end
  end

  % END output to file
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  console(log_en,'finished');
end
% END main program (control)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN function definitions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function console(enabled,string)
  if enabled
    fprintf('%.3fs:\t%s\n',toc,string);
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function symbolVectors = generate_symbolvectors(Nsym, Nsub, M, qamMode)
  switch lower(qamMode)
    case 'qam'
      Msq = sqrt(M);
      ampVals = -Msq+1:2:Msq-1;
      re_ind = randi(Msq, 1, Nsym*Nsub);
      im_ind = randi(Msq, 1, Nsym*Nsub);
      re = ampVals(re_ind);
      im = ampVals(im_ind);
      norm = sqrt(3/(2*(M-1)));
      inputSymbols = (re + 1i.*im).*norm;
    case 'pam'
      ampVals = -M+1:2:M-1;
      index = randi(M, 1, Nsym*Nsub);
      input = ampVals(index);
      norm = sqrt(3/(M^2-1));
      inputSymbols = input.*norm;
  end
  symbolVectors = reshape(inputSymbols, Nsub, length(inputSymbols)/Nsub);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function C = calculate_aic_matrix(fcz, lcz, Nsub, Nc, ovsamp)
    if lcz > fcz
        [K, L] = meshgrid(0:Nsub-1, fcz*ovsamp:lcz*ovsamp);
    else
        [K, L] = meshgrid(0:Nsub-1, [0:lcz*ovsamp fcz*ovsamp:Nsub*ovsamp-1]);
    end
    L = L/ovsamp;
    B = (K-L == 0);
    P1 = B+~B.*((1/Nsub)*((exp(2j*pi*(K-L))-1)./(~B.*exp(2j*pi/Nsub*(K-L))-1)));
  
    if lcz > fcz
        P2 = P1(:,mod((fcz-Nc:lcz+Nc),Nsub)+1);
    else
        P2 = P1(:,[0:lcz+Nc fcz-Nc:Nsub-1]+1);
    end

    m = 0;
    if m == 0
        W = pinv(P2);
    else
        PH = P2';
        S = PH*P2;
        W = inv(m*eye(length(S))+S)*PH;
    end
    C = -W*P1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function signal = calculate_OFDM_signal(s, Nsym, ecz, fcz, lcz, Nsub, Nc, opt, C)
  % initialize resulting signal with zeros
  signal = zeros(Nsub*Nsym,1);
  
  % definition of index ranges
  if lcz > fcz
      if opt.activateAIC
        cancellationIndex = mod(fcz-Nc:lcz+Nc,Nsub)+1;
      else
        cancellationIndex = mod(fcz:lcz,Nsub)+1;
      end
  else
      if opt.activateAIC
        cancellationIndex = [0:lcz+Nc fcz-Nc:Nsub-1]+1;
      else
        cancellationIndex = [0:lcz fcz:Nsub-1]+1;
      end
  end
  if opt.disableNotchCarriersNeg
    cancellationIndexNeg = (mod(-lcz,Nsub):mod(-fcz,Nsub))+1;
    cancellationIndex = [cancellationIndex  cancellationIndexNeg];
  end
  if opt.disableEdgeCarriers
    edgeCarrierIndex = mod([1:ecz Nsub-ecz:Nsub]-Nsub/2,Nsub)+1;
  end
  signalIndex = setdiff(1:Nsub, cancellationIndex);
  for k = 1:Nsym
    % get row of samples for OFDM 
    s_k = s(:, k); 
    % set DC tone to zero
    if opt.disableZeroCarrier
      s_k(1)=0;
    end
    % turn of tones at the edges
    if opt.disableEdgeCarriers
      s_k(edgeCarrierIndex) = 0;
    end
    if opt.disableNotchCarriers
      % turn of tones in the cancellation range
      s_k(cancellationIndex) = 0;
    end
    if opt.activateAIC
      % set cancellation coeffs in cancellation range
      s_k(cancellationIndex) = C*s_k;
      if opt.onlyAIC
        % deactivate the signal for better view of cancellation signal
        s_k(signalIndex) = 0;   
      end
    end
    
    % modulate OFDM symbols
    x = ifft(s_k).*sqrt(Nsub);
    % append modulated singal x to resulting signal
    signal((1:Nsub)+(k-1)*Nsub) = x;
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = generate_fir_filter(Nsub, fcz, lcz, fir_params, mode)
  % fir_params = [passband_db stopband_db nb_taps transition_width neg_en]
  passband_db = fir_params(1);
  stopband_db = fir_params(2);
  N           = fir_params(3);
  trans_width = fir_params(4);
  neg_en      = fir_params(5);
  
  if fcz > lcz
      fcz = fcz - 64;
  end
  
  gap_width = lcz-fcz;
  offset = (1+(fcz+lcz)/Nsub)/2;
  
  omega_s = (Nsub-gap_width)/Nsub;
  omega_p = (Nsub-gap_width-2*trans_width)/Nsub;
  
  omega_c = (omega_s+omega_p)/2;
  delta_f = (omega_s-omega_p)/2;
  
  switch lower(mode)
    case 'firrc'
      lowpass = generate_raised_cosine(N, omega_c, delta_f);
      h = freq_shift(lowpass, offset);
      if neg_en
        h_neg = freq_shift(lowpass, -offset);
        h = conv(h, h_neg);
      end
    case {'firpm', 'firls'}
      if    (N == 0 && stopband_db ~= 0 && passband_db ~= 0)
        % fixed passband dampening, calculate number of taps N
        delta_p = 10^(passband_db/20)-1;
        delta_s = 10^(-stopband_db/20);
        % approximation formula of bellanger (modified for only odd numbers)
        N = 2*ceil((1/3)*log10(1/(10*delta_p*delta_s))*(1/delta_f))+1;
      elseif(N ~= 0 && stopband_db == 0 && passband_db ~= 0)
        % fixed number of taps N, calculate passband dampening
        delta_p = 10^(passband_db/20)-1;
        delta_s = (1/delta_p)*10^(-(1.5*N*delta_f)-1);
      else
        delta_p = 1;
        delta_s = 1;
      end
      
      if ~neg_en
        % create complex band stop filter with real lowpass and freq. shift
        f = [0 omega_p omega_s 1];
        a = [1 1 0 0];
        w = [delta_s delta_p];
        lowpass = generate_optimized_fir(N, f, a, w, mode);
        h = freq_shift(lowpass, offset);
      else
        % create real band stop filter
        if offset > 1
            offset = 1 - rem(offset, 1);
        end
        p1 = omega_p-2*(1-offset);
        s1 = omega_s-2*(1-offset);
        s2 = s1+2*(1-omega_s);
        p2 = p1+2*(1-omega_p);
        f = [0 p1 s1 s2 p2 1]
        a = [1 1 0 0 1 1];
        w = [delta_s delta_p delta_s];
        h = generate_optimized_fir(N, f, a, w, mode);
      end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = generate_optimized_fir(N, f, a, w, mode)
  switch lower(mode)
    case 'firpm'
      h = firpm(N-1, f, a, w);
    case 'firls'
      h = firls(N-1, f, a ,w);
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = generate_raised_cosine(N, omega_c, delta_f)
  n = (-N/2:N/2-1)+0.5;
  b = abs(n) ~= 1/(2*delta_f);
  h = (b.*(cos(pi*delta_f*n)./(b.*(1-(2*delta_f*n).^2)+~b.*1))+~b*pi/4);
  h = omega_c*sinc(omega_c*n).*h;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function signal = freq_shift(signal, f)
  n = 0:length(signal)-1;
  signal = signal.*exp(2j*pi*f.*n);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function signal = apply_iq_imb(signal, phi_imb_deg, a_imb, dc_imb)
  phi_imb = pi*phi_imb_deg/180;
  
  s_re = real(signal);
  s_im = imag(signal);
  
  s_re = s_re   - a_imb*sin(phi_imb)*s_im;
  s_im = dc_imb + a_imb*cos(phi_imb)*s_im;
  
  signal = s_re + 1j*s_im;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [res, f_axis] = calculate_spectrum(signal, Nsub, Nup, dB_en, scaling)
  res_len = Nup*Nsub;
  N_nyq = res_len*ceil(length(signal)/res_len);
  if nargin == 4
    %scaling = 1/sqrt(N_nyq);
    scaling = 1/sqrt(length(signal));
  end
  y = fft(signal, N_nyq).*scaling;
  y_magn_sq = abs(fftshift(y)).^2;
  y_len = length(y_magn_sq);
  
  Nsum = y_len/res_len;
  res = zeros(1, res_len);
  for index = 1:res_len
    n_start = (index-1)*Nsum+1;
    n_end = index*Nsum;
    range = n_start:n_end;
    range_next = mod(range,y_len)+1;
    res(index) = (sum(y_magn_sq(range))+sum(y_magn_sq(range_next))) / (2*Nsum);
  end
  f_axis = ((-res_len/2:res_len/2-1)+0.5)/Nup;
  if dB_en
    res = 10*log10(res);
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function draw_points(x_axis, min_y, color)
  plot(x_axis, ones(1,length(x_axis)).*min_y,'Color',color,...
      'LineStyle','none','Marker', '.','MarkerSize', 8);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = to_Nsub_half(input, Nsub)
  result = mod(input+Nsub/2, Nsub)-Nsub/2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = get_name(lang, key)
  if strcmpi(lang, 'de')
    switch lower(key)
      case 'xlabel'
        result = 'Subcarrierindex';
      case 'firrc'
        result = 'Raised Cosine FIR ';      
      case 'firls'
        result = 'Least Squares FIR ';      
      case 'firpm'
        result = 'Equiripple FIR ';
    end
  else
    switch lower(key)
      case 'xlabel'
        result = 'subcarrierindex';
      case 'firrc'
        result = 'raised cosine FIR ';      
      case 'firls'
        result = 'least squares FIR ';      
      case 'firpm'
        result = 'equiripple FIR ';
    end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = get_legend(lang, opt, plot_2_en)
  if strcmpi(lang, 'de') 
    result = get_legend_de(opt, plot_2_en);
  else
    result = get_legend_en(opt, plot_2_en);
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function string = get_legend_de(opt, plot_2_en)
  if plot_2_en
    if opt.activateAIC
      string = {'Spektrum (AIC+Nulling) ', 'Cancellation Signal '};
    elseif opt.applyFIR
      string = {'Spektrum (FIR+Nulling) ', 'FIR Filter '};
    elseif opt.disableNotchCarriers
      string = 'Spektrum (Nulling) ';
    else
      string = 'Spektrum ';
    end
  else
    if opt.disableNotchCarriers
      string = 'Spektrum (Nulling) ';
    else
      string = 'Spektrum ';
    end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function string = get_legend_en(opt, plot_2_en)
  if plot_2_en
    if opt.activateAIC
      string = {'spectrum (AIC+nulling) ', 'cancellation signal '};
    elseif opt.applyFIR
      string = {'spectrum (FIR+nulling) ', 'FIR filter '};
    elseif opt.disableNotchCarriers
      string = 'spectrum (nulling) ';
    else
      string = 'spectrum ';
    end
  else
    if opt.disableNotchCarriers
      string = 'spectrum (nulling) ';
    else
      string = 'spectrum ';
    end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% END function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%