% script basic_noise_analysis
%
% This script evaluates 50 ohm term and receive only data.
% It loads all receive channels and then analyzes these data.
%
% For P-3:
% 1. Collect data with ? waveforms in whatever noise configuration
%    you want to measure.
% 2. Characterization should be done for 50 ohm and receive only at
%    least.
% 3. If transmitting, time gate should be large enough to include
%    noise-only data.
%
% Author: John Paden

physical_constants;
close all;
tstart = tic;

% =======================================================================
% User Settings
% =======================================================================

fs = 1e9/9;

% .plot_en = flag to enable plots
param.plot_en = true;

% .caxis = Color axis limits (leave empty first time since this causes it to use the
% defaults).
param.caxis = [];
%param.caxis = [50 120];

% .ylim = Leave empty the first time (it just uses the default limits then)
param.ylim = [];
%param.ylim = [1 500];

% .xlim = Leave empty the first time (it just uses the default limits then)
%   These limits are in range lines post presumming
param.xlim = [];
%param.xlim = [200 450];

% .ref_adc = Reference receive channel (surface location determined from this
%   channel and all phase measurements made relative to it)
param.ref_adc = 5;
% .rlines = Range lines to process
%   These are range lines post presumming
param.rlines = [1 inf];
% .rbins = Range bins to search for surface in
param.rbins = [1 inf];
% .noise_rbins = Range bins to use for noise power calculation (THIS OFTEN NEEDS TO BE SET)
param.noise_rbins = 4000:7000;

% .snr_threshold = SNR threshold in dB (range lines exceeding this
%   SNR are included in the estimate)
param.snr_threshold = 25;

param.wf = 2;
param.pc_param.f0 = 180e6;
param.pc_param.f1 = 210e6;
param.pc_param.Tpd = 1e-6;
param.pc_param.tukey = 0.2;
param.pc_param.rx_gain = 1;

param.radar_name = 'mcords';
if strcmpi(param.radar_name,'mcords')
  param.radar_num = 3;
  param.adcs = [1 2 3 4 5 6];
  param.seg = 1;
  
  % Parameters to locate specific file of interest
  %    (THIS NEEDS TO BE SET EVERYTIME)
  param.data_file_num = 3;
%   param.base_path = '/home/polargrid/mcords/2010_Antarctica_DC8/20101106/';
  param.base_path = '/mnt/mcords_new1/';
elseif strcmpi(param.radar_name,'mcords2')
  % .adcs = the receive channels to use
  param.adcs = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
  
  param.base_path = 'D:\mcords2\2011_Greenland_P3\20110314\';
  param.acquisition_num = [];
  param.seg = 13;
  param.file_num = 1;
end

% .gps_fn = Optional GPS file name (leave empty to disable)
param.gps_fn = '/home/polargrid/csarp_support/gps/2010_Antarctica_DC8_GPS/20101026_ALL_pos.mat';
param.gps_fn = '';

% .presums = Number of presums (coherent averaging) to do
param.presums = 50;

% .noise_removal_en = mcords noise removal (i.e. should be false for
%    mcrds, mcords2, etc)
param.noise_removal_en = false;

adc_bits = 14;
Vpp_scale = 2;
adc_SNR_dB = 70;
rx_gain = 10^((72-20)/20);
noise_figure = 10^(1.6/10); % Do not include receiver losses

rline = 1;

% =======================================================================
% =======================================================================
% Automated Section
% =======================================================================
% =======================================================================

% =======================================================================
% Load data
% =======================================================================
fprintf('========================================================\n');
fprintf('Loading data (%.1f sec)\n', toc(tstart));
% Load the data (disable if you have already loaded)
clear data;
clear num_rec;
if strcmpi(param.radar_name,'mcords')
  for adc_idx = 1:length(param.adcs)
    adc = param.adcs(adc_idx);

    % May need to adjust base_path for non-standard directory structures
%     base_path = fullfile(param.base_path, sprintf('chan%d',adc), ...
%       sprintf('seg_%d',param.seg));
    base_path = fullfile(param.base_path, sprintf('mcords_new%d',adc));
    file_prefix = sprintf('mcords');
    file_midfix = sprintf('r%d-%d.',param.radar_num,adc);
    file_suffix = sprintf('.%04d.dat',param.data_file_num);
    fprintf('  Path: %s\n', base_path);
    fprintf('  Match: mcords*%s*%s\n', file_midfix, file_suffix);
    fn = get_filename(base_path,'mcords',file_midfix,file_suffix);
    if isempty(fn)
      fprintf('  Could not find any files which match\n');
      return;
    end
    [hdr,data] = basic_load_mcords(fn, struct('clk',1e9/9,'first_byte',2^26));

    fn = get_filename(base_path, file_prefix, '', file_suffix);
    if isempty(fn)
      fprintf('  Could not find any files which match\n');
      return;
    end
    fprintf('  Loading file %s\n', fn);
    [hdr,data_tmp] = basic_load_mcords(fn,struct('clk',fs));
%     data(:,:,adc_idx) = data_tmp{param.wf};
    data{adc_idx} = data_tmp{param.wf};
  end
  basic_remove_mcords_digital_errors;
elseif strcmpi(param.radar_name,'mcords2')
  % test1_1.dat0
  %   testA_N.datB
  %   A = acquisition number
  %   N = file number
  %   B = board number
  % Map ADCs to board numbers
  for board = 0:3
    if any(board == floor((param.adcs-1)/4))
      get_adcs = board*4 + (1:4);
      file_prefix = sprintf('mcords2_%d_',board);
      if isempty(param.acquisition_num)
        file_suffix = sprintf('%04d.bin',param.file_num);
      else
        file_suffix = sprintf('%02d_%04d.bin',param.acquisition_num,param.file_num);
      end
      base_path = fullfile(param.base_path, sprintf('board%d',board), ...
        sprintf('seg_%02d',param.seg));
      fprintf('  Path: %s\n', base_path);
      fprintf('  Match: %s*%s\n', file_prefix, file_suffix);
      fn = get_filename(base_path, file_prefix, '', file_suffix);
      if isempty(fn)
        fprintf('  Could not find any files which match\n');
        return;
      end
      fprintf('  Loading file %s\n', fn);
      % Fix get_filenames     'The filename, directory name, or volume label syntax is incorrect.'
      [hdr,data_tmp] = basic_load_mcords2(fn,struct('clk',fs));
      for get_adc_idx = 1:length(get_adcs)
        for adc_idx = 1:length(param.adcs)
          if param.adcs(adc_idx) == get_adcs(get_adc_idx)
            if ~exist('num_rec','var')
              % Since each file may have slightly different numbers of
              % records we do this
              num_rec = size(data_tmp{param.wf},2) - 1;
            end
            data(:,:,adc_idx) = data_tmp{param.wf}(:,1:num_rec,get_adc_idx);
          end
        end
      end
    end
  end
end
[fn_path fn_name] = fileparts(fn);
if ~isfield(param,'seg') || isempty(param.seg)
  param.seg = -1;
end
clear data_tmp;

% =======================================================================
% Convert from quantization to voltage @ ADC
data = data ...
  * Vpp_scale/2^adc_bits ...
  * 2^hdr.wfs(param.wf).bit_shifts/hdr.wfs(param.wf).presums;

% =====================================================================
% Noise power
% =====================================================================

% Calculate noise power as Vrms
fprintf('Expected ADC noise floor @ ADC %.1f dBm\n', lp((Vpp_scale/2/sqrt(2))^2/50)+30 - adc_SNR_dB );
fprintf('Expected Rx noise floor @ ADC %.1f dBm\n', lp(BoltzmannConst*290*30e6*noise_figure*rx_gain^2)+30);
fprintf('Noise power is in dBm:\n')
for adc_idx = 1:size(data,3)
  if adc_idx > 1
    fprintf('\t');
  end
  fprintf('%.1f', ...
    lp(mean(mean(abs(data(param.noise_rbins,:,adc_idx)).^2*hdr.wfs(param.wf).presums/50)))+30 );
end
fprintf('\n');

% =====================================================================
% Quantization analysis
% =====================================================================
if 0
  for adc_idx = 1:size(data,3)
    data = data - mean(data(:,rline));
    
    figure(adc_idx); clf;
    plot(data(:,rline,adc_idx),'.');
    grid on;
    xlabel('Range bin');
    ylabel('Quantization level');
    
    figure(20+adc_idx);
    imagesc(lp(data(:,:,adc_idx),2));
    colorbar;
    
    % Plot estimated pdf and approximate a gaussian to it
    figure(40+adc_idx);
    ROI = data(param.noise_rbins,100:200,adc_idx);
    ROI = ROI(:);
    [n,x] = hist(ROI,64);    
    bar(x,n);
    mean_x = mean(ROI);
    var_x = var(ROI);
    hold on;
    plot(x, numel(ROI)*(x(2)-x(1)) * 1/sqrt(2*pi*var_x) * exp(-(x - mean_x).^2 / (2*var_x)),'r');
    hold off;
  end
  for adc_idx = 1:size(data,3)
    set(adc_idx,'WindowStyle','docked','NumberTitle','off','Name',sprintf('Q%d',adc_idx));
  end  
  for adc_idx = 1:size(data,3)
    set(20+adc_idx,'WindowStyle','docked','NumberTitle','off','Name',sprintf('E%d',adc_idx));
  end  
  for adc_idx = 1:size(data,3)
    set(40+adc_idx,'WindowStyle','docked','NumberTitle','off','Name',sprintf('P%d',adc_idx));
  end  
end

% =====================================================================
% Power Spectrum
% =====================================================================
if 1
  for adc_idx = 1:size(data,3)
    adc = param.adcs(adc_idx);
    
    if ~isfield(param,'noise_rbins') || isempty(param.noise_rbins)
      fir_data = fir_dec(data(:,:,adc_idx),param.presums);
    else
      fir_data = fir_dec(data(param.noise_rbins,:,adc_idx),param.presums);
    end
    
    pc_param.time = hdr.wfs(param.wf).t0 + (0:size(fir_data,1)-1)/fs;
    dt = pc_param.time(2) - pc_param.time(1);
    Nt = length(pc_param.time);
    df = 1/(Nt*dt);
    freq = fs + (0:df:(Nt-1)*df).';
        
    figure(100+adc); clf;
    set(100+adc,'WindowStyle','docked','NumberTitle','off','Name',sprintf('FFT%d',adc_idx));
    imagesc([], freq/1e6, lp(fft(fir_data)) + 30 + 10*log10(2^2/50/size(fir_data,1)) )
    title(sprintf('Freq-space adc%d ave%d seg%d/%s', adc, param.presums, param.seg, fn_name),'Interpreter','none');
    xlabel('Range line');
    ylabel('Frequency (MHz)');
    ylim(fs/1e6*[1.5 2]);
    h = colorbar;
    set(get(h,'YLabel'),'String','Relative power (dB)');
    
    figure(120+adc); clf;
    set(120+adc,'WindowStyle','docked','NumberTitle','off','Name',sprintf('M%d',adc_idx));
    plot(freq/1e6, lp(mean(abs(fft(fir_data)).^2*2^2 / 50,2)/size(fir_data,1)) + 30)
    title(sprintf('MeanFFT adc%d ave%d seg%d/%s', adc, param.presums, param.seg, fn_name),'Interpreter','none');
    ylabel('Relative power (dB)');
    xlabel('Frequency (MHz)');
    xlim(fs/1e6*[1.5 2]);
    grid on;
  end
end









