% spectrumAnalysis
%

clear;clc;
%% Setup
addpath(genpath('../../gss/navsol'));
datadir = './turkey_captures/';
filename = 'capture_20231123_005';
% Interval of data to load (sec)
Tp = 20/750;
% Seek time into data (sec)
tSeek = 53;
% Size of FFT used in power spectrum estimation
Nfft = 2^12;        
% Indicates whether xml file is available
xmlAvailable = true;
% Indicate whether to choose power spectrum (p) or spectrogram (s)
displayFlag = 's';

%% Load data
if(xmlAvailable)
[t0, Fsr, Fcr] = parseNIxml(datadir, filename);
t0 = t0 + seconds(tSeek);
% t0 is expressed in UTC; convert to GPS time
[t0_gps_week, t0_gps_sec] = datetime2gps(utc2gps(t0));
else
  % Substitute capture details if xml not available
  t0_gps_week = 0;
  t0_gps_sec = 0;
  Fsr = 50e6;
  Fcr = 12.050e9;
end
[yVec,tyVec] = loadNIdata(datadir,filename,Tp,Fsr,tSeek);
Ny = length(yVec);

%% Filter
if(1)
  f3dB = 28.2e6;
  orderLowpass = 8;
  d = fdesign.lowpass('N,F3dB,Ap',orderLowpass,f3dB,0.5,Fsr);
  Hd = design(d, 'cheby1');
  yVec = filter(Hd,yVec);
end

%% Frequency shift 
% The center frequency of the ith Starlink band is given by Fci = 10.7 + F/2 +
% 0.25(i-1/2), in GHz.  Frequency shift the data to account for the difference
% between Fci and Fcr.
Fdelta = 250e6;
Fs = 240e6;
N = 1024;
F = Fs/N;
ii = round((Fcr/1e9 - 10.7 - F/2/1e9)/0.25 + 0.5);
Fcii = (10.7e9 + F/2 + Fdelta*(ii - 0.5));
fShift0 = Fcr - Fcii;
tyVec = [0:Ny-1]'/Fsr;
yVec = yVec.*exp(j*2*pi*fShift0*tyVec);

%% Compute power time history
figure(1);clf;
Nsamps_per_block = 0.0002*Fsr;
Nblocks = floor(Ny/Nsamps_per_block);
M = zeros(Nsamps_per_block,Nblocks);
yVec = yVec(1:Nsamps_per_block*Nblocks);
M(:) = yVec.*conj(yVec);
PVec = mean(M,1);
tVec = [0:Nsamps_per_block:Nsamps_per_block*(Nblocks-1)]*(1/Fsr);
area(tVec,10*log10(PVec),42);
xlabel('Time since t_0 (sec)');
ylabel('Power (dB)');
title('Power time history');
grid on;

if(displayFlag == 'p')
  %% Compute power spectrum estimate
  [Syy,fVec] = pwelch(yVec,kaiser(Nfft,3),Nfft/2,Nfft,Fsr,'psd','centered');

  %% Plot results
  figure(2);clf;
  plot(fVec/1e6,10*log10(Syy));
  grid on;
  xlabel('Frequency (MHz)');
  ylabel('Power density (dB/Hz)');
  title('Power spectral density estimate');
  ylim([-35 -10]);
else
  %% Compute spectrogram
  figure(2);clf;
  spectrogram(yVec,kaiser(Nfft,0.5),Nfft/2,Nfft,Fsr,'psd','centered','yaxis', ...
              'minthreshold', -30);
  title('Spectrogram');
  grid on;
end

commandwindow


  
