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 = 0;
  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 bssc = data(:,10);
  55 bsic = data(:,11);
  56 bsid = data(:,10)*1000 + data(:,11);
  57 % Received signal level
  58 rxl = data(:,12);
  59 % MCC/MNC numbers for the network
  60 mcc = data(1,6);
  61 mnc = data(1,7);
  62
  63 %% Getting all distinct CellID numbers
  64 bsidlist = unique(bsid);
  65 %% Sometimes there's unknown cellid logged as '0' - omit
  66 if (bsidlist(1) == 0)
  67     bsidlist = bsidlist(2:end);
  68 end % if
  69 disp(['Number of cells in log: ' num2str(length(bsidlist))])
  70 disp(['MCC,MNC: ' num2str(mcc) ',' num2str(mnc)]);
  71
  72 if (cellcheck < 1) or (cellcheck > length(bsidlist))
  73         disp('Invalid cell number - no such cell to check')
  74 end % if
  75
  76
  77 % plot complete logged GPS track
  78 if toplot == 1
  79         hold off
  80         figure(1)
  81         plot(lon,lat,'g.')
  82         axis equal
  83         hold on
  84 end % if
  85
  86 out = [];
  87
  88 %% Do fitting
  89 %%%%% In loop to allow for later extending to batch procession
  90 %%%%% Currently only one at a time: need of visual inspection of data for each cell
  91 for k = cellcheck:cellcheck
  92
  93         %% Current CellID
  94         id = bsidlist(k);
  95         bssc_k = floor(id/1000);
  96         bsic_k = id - bssc_k*1000;
  97         disp(['MCC,MNC,Combined ID,BSIC,BSSC: ' num2str(mcc) ',' num2str(mnc) ',' num2str(id) ',' num2str(bsic_k) ',' num2str(bssc_k) ])
  98
  99         %% Number of logged points
 100         nlog = length(bsid(bsid == id));
 101         disp(['Logged points: ' num2str(nlog)]);
 102         if nlog == 1
 103                 disp('Only one point in data, cannot calculate anything');
 104                 exit;
 105         end
 106
 107         %% Select logged parameters for just this cell
 108         clat =  lat(bsid == id);
 109         clon =  lon(bsid == id);
 110         crxl =  rxl(bsid == id);
 111         cbs = bsic(bsid == id);
 112         basestations = unique(cbs);
 113         disp(["Base stations: " num2str(length(basestations))])
 114
 115         %% Turn RxL (dB) into signal scale
 116         csig =  10.^(crxl/10);
 117
 118         % Calculate all signal compared to maximum recorded
 119         crxlnorm = max(crxl)-crxl;
 120         % Distance scaling power: d^(-dsq) : dsq=2 is free space propagation
 121         dsc = 2;
 122         %%%% scale factor only for plotting, that should be handled more cleverly...
 123         dscale = 0.00001;
 124
 125         %% Approximate distance scale:
 126         % RxL (dB) -> Signal (relative) -> Distance (with assumed scaling)
 127         ds = 1./(10.^(crxlnorm/10)).^(1/dsc);
 128         dist = dscale*1./ds;
 129
 130         %% Plot data recorded for given cell
 131         % Location helper:
 132         %   - circles around logged positions
 133         %   - circle radius is inversely proportional to signal strength
 134         %   ---> tower 'close': small circle; tower 'far': big circle
 135         if toplot == 1
 136                 figure(1)
 137                 plot(clon,clat,'b.')
 138                 for l = 1:nlog
 139                         c1 = dist(l)*cos(linspace(0,2*pi))+clon(l);
 140                         c2 = dist(l)*sin(linspace(0,2*pi))+clat(l);
 141                         plot(c1,c2)
 142                 end % for
 143                 axis equal
 144         end % if
 145
 146
 147         %% For fitting, guess the tower position close to the max recorded signal
 148         %% an add random offset
 149         [maxsig maxsignum] = max(crxl);
 150         offsetm = 0.003;
 151         fitofflon = offsetm*rand-offsetm/2;
 152         fitofflat = offsetm*rand-offsetm/2;
 153         b_in = [clon(maxsignum)+fitofflon clat(maxsignum)+fitofflat 0.01];
 154
 155         %% use the RxL for fitting...
 156         [yhat, b_est] = leasqr([clon clat],crxl,b_in,"gsmsurf",0.00001,3000);
 157         disp(['Fitted tower coordinates (Lat/Lon):  ' num2str(b_est(2)) ',' num2str(b_est(1))]);
 158
 159         % Residual sum of squares (for fit quality purposes: smaller the better!
 160         RSSin = sum((gsmsurf([clon clat],b_in)-crxl).^2);
 161         RSSest = sum((gsmsurf([clon clat],b_est)-crxl).^2);
 162         disp(['RSS: ' num2str(RSSin) ' -> ' num2str(RSSest)])
 163
 164         %%% Put fitted tower position on map: big red X
 165         if toplot == 1
 166                 plot(b_est(1),b_est(2),'rx','MarkerSize',32);
 167                 xlabel(['CellID:' num2str(id)  '/ RSS: ' num2str(RSSest) ' / Pos: ' num2str([b_est(1) b_est(2)])])
 168                 axis equal;
 169                 hold off
 170         end % if
 171
 172         %% To assess fit quality:
 173         if toplot == 1
 174                 figure(2)
 175                 lon0 = b_est(1);
 176                 lat0 = b_est(2);
 177                 r = sqrt((clon-lon0).^2 + (clat-lat0).^2);
 178 %               yhat = yhat./log(10)*10;
 179                 plot(r,crxl,'@',r,yhat,'x')
 180 %               plot(r,lcsig,'@',r,yhat,'x')
 181                 ylabel('RxL (dB)')
 182                 xlabel(['Distance from fitted tower position / CellID:' num2str(id) ])
 183         end % if
 184
 185         disp(['GSMPOS: ' num2str(mcc) ',' num2str(mnc) ',' num2str(id) ',' num2str(b_est(2)) ',' num2str(b_est(1)) ])
 186         disp(['OpenLayers: ' num2str(b_est(2)) '|' num2str(b_est(1)) '|cell ' num2str(id) '|mcc ' num2str(mcc) ' mnc ' num2str(mnc) '|icon_blue.png|24,24|0,-24'])
 187 end % for
 188 end % function gsmfit
 189
 190
 191 %%%%%%
 192 %%% Additional sub-routines
 193 %%%%%%
 194
 195 function yhat = gsmsurf(X,b)
 196 %% Signal propagation surface according to propagation assumption
 197
 198         %% Assumed distance scaling: d^(-dsc): dsc=2 is free-space propagation
 199         dsc = 2;
 200         %%%%% Separating input parameters into variables
 201         % proposed cell position
 202         lon0 = b(1);
 203         lat0 = b(2);
 204         % tower output scale factor - "signal at tower" -  purely for fitting
 205         s0 = b(3);
 206         %%%%%
 207
 208         % Distance from checked centre
 209         r = sqrt((X(:,1)-lon0).^2 + (X(:,2)-lat0).^2);
 210         % Estimated 'signal'
 211         yhat = 10*log(s0*r.^(-dsc))/log(10);
 212 end
 213
 214 %%%%%%%%% END
 215
 216 gsmfit