gsmpos/gsmfit.m

   1 %%%%%%%%%%%%%%%%%%
   2 %  (C) 2008 Gergely Imreh  <imrehg@gmail.com>
   3 %
   4 %  GPLv3 or later
   5 %%%%%%%%%%%%%%%%%%
   6
   7 %%%%%%%%%%%%%%%%%%
   8 % Crude Octave3 code to attempt to find GSM tower positions
   9 % using log from Cell Locator:
  10 % http://repo.or.cz/w/CellLocator.git
  11 % http://projects.openmoko.org/projects/rocinante/
  12 %
  13 % Very rough, trying to fit free-space propagation model to
  14 % recorded signal. Also visualizing log data.
  15 %%%%%%%%%%%%%%%%%%
  16
  17 function gsmfit()
  18
  19 %%%%%%%%%%%%
  20 % Parameters to modify
  21 %%%%%%%%%%%%
  22 % GSM/GPS Logfile
  23 logfile = 'log_cells.csv';
  24 % Clean up before use! Remove headers, and possible invalid lines
  25
  26 % Which GSM cell to check
  27 cellcheck = 1;
  28 % it is the number in the order of the cells listed in logfile
  29 % cellcheck = [1..number-of-different-cellid]
  30
  31 % Whether to plot: 0/1
  32 toplot = 1;
  33
  34 %%%%End parameters
  35
  36 disp('Parameters ==============')
  37 disp(['Logfile: ',logfile])
  38 disp(['Cell to check : ' num2str(cellcheck)])
  39
  40 close all;
  41 format none
  42
  43 %% Load logfile
  44 data = load(logfile);
  45
  46
  47 %% Parsing log into separate variables
  48 % GPS position
  49 lon = data(:,3);
  50 lat = data(:,4);
  51 % CellID
  52 cid = data(:,9);
  53 % Base station id
  54 bs = data(:,11);
  55 % Received signal level
  56 rxl = data(:,12);
  57 % MCC/MNC numbers for the network
  58 mcc = data(1,6);
  59 mnc = data(1,7);
  60
  61 %% Getting all distinct CellID numbers
  62 cidlist = unique(cid);
  63 % Sometimes there's unknown cellid logged as '0' - omit
  64 if (cidlist(1) == 0)
  65     cidlist = cidlist(2:end);
  66 end % if
  67 disp(['Number of cells in log: ' num2str(length(cidlist))])
  68 disp(['MCC,MNC: ' num2str(mcc) ',' num2str(mnc)]);
  69
  70 if (cellcheck < 1) or (cellcheck > length(cidlist))
  71         disp('Invalid cell number - no such cell to check')
  72 end % if
  73
  74
  75 % plot complete logged GPS track
  76 if toplot == 1
  77         hold off
  78         figure(1)
  79         plot(lon,lat,'g.')
  80         axis equal
  81         hold on
  82 end % if
  83
  84 out = [];
  85
  86 %% Do fitting
  87 %%%%% In loop to allow for later extending to batch procession
  88 %%%%% Currently only one at a time: need of visual inspection of data for each cell
  89 for k = cellcheck:cellcheck
  90
  91         %% Current CellID
  92         id = cidlist(k);
  93         disp(['MCC,MNC,CellID: ' num2str(mcc) ',' num2str(mnc) ',' num2str(id)])
  94
  95         %% Number of logged points
  96         nlog = length(cid(cid == id));
  97         disp(['Logged points: ' num2str(nlog)])
  98
  99         %% Select logged parameters for just this cell
 100         clat =  lat(cid == id);
 101         clon =  lon(cid == id);
 102         crxl =  rxl(cid == id);
 103         cbs = bs(cid==id);
 104         basestations = unique(cbs);
 105         disp(["Base stations: " num2str(basestations')])
 106
 107         %% Turn RxL (dB) into signal scale
 108         csig =  10.^(crxl/10);
 109
 110         % Calculate all signal compared to maximum recorded
 111         crxlnorm = max(crxl)-crxl;
 112         % Distance scaling power: d^(-dsq) : dsq=2 is free space propagation
 113         dsc = 2;
 114         %%%% scale factor only for plotting, that should be handled more cleverly...
 115         dscale = 0.00001;
 116
 117         %% Approximate distance scale:
 118         % RxL (dB) -> Signal (relative) -> Distance (with assumed scaling)
 119         ds = 1./(10.^(crxlnorm/10)).^(1/dsc);
 120         dist = dscale*1./ds;
 121
 122         %% Plot data recorded for given cell
 123         % Location helper:
 124         %   - circles around logged positions
 125         %   - circle radius is inversely proportional to signal strength
 126         %   ---> tower 'close': small circle; tower 'far': big circle
 127         if toplot == 1
 128                 figure(1)
 129                 plot(clon,clat,'b.')
 130                 for l = 1:nlog
 131                         c1 = dist(l)*cos(linspace(0,2*pi))+clon(l);
 132                         c2 = dist(l)*sin(linspace(0,2*pi))+clat(l);
 133                         plot(c1,c2)
 134                 end % for
 135                 axis equal
 136         end % if
 137
 138
 139         %% For fitting, guess the tower position close to the max recorded signal
 140         %% an add random offset
 141         [maxsig maxsignum] = max(crxl);
 142         offsetm = 0.003;
 143         fitofflon = offsetm*rand-offsetm/2;
 144         fitofflat = offsetm*rand-offsetm/2;
 145         b_in = [clon(maxsignum)+fitofflon clat(maxsignum)+fitofflat 0.01];
 146
 147         %% use the RxL for fitting...
 148         [yhat, b_est] = leasqr([clon clat],crxl,b_in,"gsmsurf",0.00001,3000);
 149         disp(['Fitted tower coordinates (Lat/Lon):  ' num2str(b_est(2)) ',' num2str(b_est(1))]);
 150
 151         % Residual sum of squares (for fit quality purposes: smaller the better!
 152         RSSin = sum((gsmsurf([clon clat],b_in)-crxl).^2);
 153         RSSest = sum((gsmsurf([clon clat],b_est)-crxl).^2);
 154         disp(['RSS: ' num2str(RSSin) ' -> ' num2str(RSSest)])
 155
 156         %%% Put fitted tower position on map: big red X
 157         if toplot == 1
 158                 plot(b_est(1),b_est(2),'rx','MarkerSize',32);
 159                 xlabel(['CellID:' num2str(id)  '/ RSS: ' num2str(RSSest) ' / Pos: ' num2str([b_est(1) b_est(2)])])
 160                 axis equal;
 161                 hold off
 162         end % if
 163
 164         %% To assess fit quality:
 165         if toplot == 1
 166                 figure(2)
 167                 lon0 = b_est(1);
 168                 lat0 = b_est(2);
 169                 r = sqrt((clon-lon0).^2 + (clat-lat0).^2);
 170 %               yhat = yhat./log(10)*10;
 171                 plot(r,crxl,'@',r,yhat,'x')
 172 %               plot(r,lcsig,'@',r,yhat,'x')
 173                 ylabel('RxL (dB)')
 174                 xlabel(['Distance from fitted tower position / CellID:' num2str(id) ])
 175         end % if
 176
 177         disp(['GSMPOS: ' num2str(mcc) ',' num2str(mnc) ',' num2str(id) ',' num2str(b_est(1)) ',' num2str(b_est(2)) ])
 178         disp(['OpenLayers: ' num2str(b_est(1)) '|' num2str(b_est(2)) '|cell ' num2str(id) '|mcc ' num2str(mcc) ' mnc ' num2str(mnc) '|icon_blue.png|24,24|0,-24'])
 179 end % for
 180 end % function gsmfit
 181
 182
 183 %%%%%%
 184 %%% Additional sub-routines
 185 %%%%%%
 186
 187 function yhat = gsmsurf(X,b)
 188 %% Signal propagation surface according to propagation assumption
 189
 190         %% Assumed distance scaling: d^(-dsc): dsc=2 is free-space propagation
 191         dsc = 2;
 192         %%%%% Separating input parameters into variables
 193         % proposed cell position
 194         lon0 = b(1);
 195         lat0 = b(2);
 196         % tower output scale factor - "signal at tower" -  purely for fitting
 197         s0 = b(3);
 198         %%%%%
 199
 200         % Distance from checked centre
 201         r = sqrt((X(:,1)-lon0).^2 + (X(:,2)-lat0).^2);
 202         % Estimated 'signal'
 203         yhat = 10*log(s0*r.^(-dsc))/log(10);
 204 end
 205
 206 %%%%%%%%% END