%% Quick unwrap of broad stepped frequency
clear
%filename = 'ecdr_chip02_chan00@1308961140_stripped.bin' %1MHz, 150MHz
125ms dwell
%filename = 'ecdr_chip02_chan00@1309388220_stripped.bin' %1MHz, 30MHz,
600kHz step 125ms dwell
%filename = 'ecdr_chip02_chan00@1309450740_stripped.bin'
% filename = 'ecdr_chip02_chan00@1309465140_stripped.bin'
% basepath = 'uw-isis002.ee.washington.edu/isis-rx21/isis/2011-06-30/';

filename = 'ecdr_chip02_chan00@1309566900_stripped.bin';
basepath = 'uw-isis002.ee.washington.edu/isis-rx21/isis/2011-07-02/';

filepath = [basepath,filename];

bandwidth = 1e6
%time_offset = 3.7702
%time_offset = 2.3;
time_offset = 0;

sample_offset = time_offset*bandwidth;


freq_step = 500e3
freq_span = 500e6
num_steps = freq_span/freq_step;
dwell = 125e-3

nfft = 512;
fft_res = bandwidth/nfft;
no_fft_overlap = (bandwidth-freq_step)/2/fft_res;
no_fft_overlap = floor(no_fft_overlap);
%no_fft_overlap = 1;

chunk_size = 2*dwell*bandwidth;
num_avgs = floor(chunk_size/(2*nfft));



%% Create the DPSS data tapers
num_tapers = 16;
[data_window, conc] = dpss(nfft, 4 , num_tapers);

%Create the tapering matrix.
W = repmat(data_window,1,floor(num_avgs/num_tapers));
if mod(num_avgs,num_tapers)
   W = [W,W(:, 1:mod(num_avgs,num_tapers))];
end

%Bracket the averages
NBLANK = 40;
W(:,1:NBLANK) = zeros(nfft,NBLANK);
W(:,end-NBLANK+1:end) = zeros(nfft,NBLANK);

%% Get number of samples
fid = fopen(filepath, 'rb');
fseek(fid, 0, 1);
num_samples = ftell(fid)/4;
%Rewind back to sample_offset
fseek(fid, 4*sample_offset, -1);


unwrapped_Pxx = zeros(num_steps*(nfft-2*no_fft_overlap),1);

lbound = 1;
ubound = nfft-2*no_fft_overlap;
%% Main loop
for step_no = [0:num_steps-1];
   iq_raw = fread(fid, chunk_size, 'int16');

   %Have we reached the end of the file?  Quit
   if length(iq_raw) < chunk_size
       done = true;
       break
   end

   %Unpack IQ into complex samples
   iq = iq_raw(1:2:end) + 1j*iq_raw(2:2:end);

   %Reshape into an [nfft x num_avg] matrix, take FFT along row after
   %tapering with the window matrix W
   iq = reshape(iq(1:nfft*num_avgs), nfft, num_avgs);
   iq_fft = fftshift(fft(iq.*W, nfft, 1),1);

   mean_fft = mean(iq_fft .* conj(iq_fft), 2);

   %if step_no > 170
   %step_no
   %end
   unwrapped_Pxx(lbound:ubound) = unwrapped_Pxx(lbound:ubound) + ...
       mean_fft(no_fft_overlap:end-no_fft_overlap-1);



   if mod(step_no,20) == 0
       disp(sprintf('Spectra %d of %d complete', step_no, num_steps))

   end

   %Shift the bounds
   lbound = ubound + 1;% + no_fft_overlap;
   ubound = ubound + nfft - 2*no_fft_overlap;

end
fclose(fid);
Fc = 850e6;
F = 0:1:num_steps*(nfft-2*no_fft_overlap)-1;
figure; plot((F*fft_res + Fc + no_fft_overlap*fft_res)/100e4,
10*log10(unwrapped_Pxx))