function []=calibration_filter(filein)
%%%%% Calibration for CRX filters %%%
%
% Compute the phase differences respect antenna 1 induced due the filters inside
% the CRX given an calibration acquisition in .mat
%
% The output is a 3 column matrix file where:
% each row represents an antenna
% column 1 are idealized normalized amplitude values for channel I
% column 2 are idealized normalized amplitude values for channel Q
% column 3 are the phase differences respect antenna 1 in degrees
%
% Move this calibration file to ~/etc/ and make a symbolic link called
% calibration.wera pointing to this file
%
% It doesnt matter column 1 and 2 are filled with ones because Gurgel
% software re-normalize the amplitude according to the Gain of the antennas 
%
%
% Isaac Rodriguez-Padilla march/2017
% Johanna, October/2024 modified to read the data directly from dta
%clear all; close all; clc


load /home/ruser/conf/radar_header.mat

cd /home/ruser/data/abcd/

%load(filein)

% number of samples to overlap on each side
shift_pos = floor(RC.SHIFT_FRAC.*RC.dtacq_samp_chirp/2); % make sure that it is whole

% length of chirp with added overlap samples
samp_chirp_lengthen = RC.dtacq_samp_chirp + shift_pos*2;

%filein='20242780301_abcd.dta';

% open and prepare input file
fi = fopen(filein,'r','ieee-le');
fseek(fi,0,'eof');
pos = ftell(fi);

% how many possible chirps are there in the file?
%  (after skipping the first )
file_chirps = floor((pos/4)/(RC.NCHAN*RC.IQ*RC.dtacq_samp_chirp))-2;
clear pos

fseek(fi,0,'bof'); %beginning of the file
data=fread(fi,[RC.NCHAN*RC.IQ,RC.dtacq_samp_chirp],'int32=>double');  % skip 1 chirp's worth of data

%fastforward shift_pos so that rewind in loop is consistent
fseek(fi,(RC.NCHAN*RC.IQ*(shift_pos))*4,'cof');

%for nchirp=1:2
fseek(fi,-(RC.NCHAN*RC.IQ*shift_pos*2)*4,'cof');
data = fread(fi,[RC.NCHAN*RC.IQ,samp_chirp_lengthen],'int32=>double');
%end
fclose(fi);

datai=squeeze(data(1:2:end,:));
dataq=squeeze(data(2:2:end,:));
% Reference antenna 1
ch1i = datai(1,:); % channel I, antenna 1
ch1q = dataq(1,:); % channel Q, antenna 1
comp1 = ch1i+i*ch1q; % complex number for antenna 1

% Change the loop value depending on the number of antennas 
for ANT=2:RC.NANT;

    chi = datai(ANT,:); % channel I
    chq = dataq(ANT,:); % channel Q

    comp = chi+i*chq; % complex number for each antenna

    deltaphi=angle(comp./comp1); % compute the difference of angles respect ant 1
 % deltaphi=nanmedian(comp./comp1);   
   phase_dif=median(deltaphi(:)); % compute one mean value for the phase difference. We use the median because we're not dealing with a normal distribution 
  %  phase_dif=angle(deltaphi);
    phase_diff(1)=0;
    phase_diff(ANT)=rad2deg(phase_dif);

end

format long
phase_diff' % Phase differences in degrees respect ant 1

% make the calibration matrix
cal = ones(length(phase_diff),3);
cal(:,3) = phase_diff';

cd /home/wera/wera_import
% generate a .txt file
eval(['fid = fopen(''calibration.wera_',date,''',''wt'');']);
%fid = fopen('calibration.wera','wt');
fprintf(fid,'%14.8f %14.8f %14.8f\n',cal');
fclose(fid);

% Save all parameters in .mat file for archive
%save(['rx_calib_',date,'.mat',',cal']);
fileout=['rx_calib.mat'];
save(fileout,'cal');


