From e15b1240ed4553a35ac1fb8827f0a7694ab7669a Mon Sep 17 00:00:00 2001 From: Gergely Imreh Date: Thu, 25 Sep 2008 18:22:31 +0100 Subject: [PATCH] Add: finding cell tower position from log An Octave3 script finding GSM tower position based on log file collected with Cell Locator. Very crude at the moment, but even when it cannot find the right location, one can use the plot to read off a better guess (gnuplot shows plot coordinate at the cursor's position so it can be used as a map). --- gsmpos/gsmfit.README | 55 ++++++++++++++ gsmpos/gsmfit.m | 204 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 259 insertions(+) create mode 100644 gsmpos/gsmfit.README create mode 100644 gsmpos/gsmfit.m diff --git a/gsmpos/gsmfit.README b/gsmpos/gsmfit.README new file mode 100644 index 0000000..dfd2c0a --- /dev/null +++ b/gsmpos/gsmfit.README @@ -0,0 +1,55 @@ +-------- gsmfit() ---------------------------------------- + +Crude Octave3 code to attempt to find GSM tower positions +using log from Cell Locator: +http://repo.or.cz/w/CellLocator.git +http://projects.openmoko.org/projects/rocinante/ + +Very rough, trying to fit free-space propagation model to +recorded signal. Also visualizing log data. + +---------------------------------------------------------- + +To use, change the following variables in the source: + +logfile: +-------- +The path to the log to be used + +cellcheck: +---------- +Which cell should be filtered for in the logfile. This is NOT +the CellID, but [1...number-of-unique-cellid-in-log]. +Right now the analysis is not reliable enough to allow for batch +processing. That should come maybe in future editions. + +toplot: +------- +0/1: No plotting / plotting of logged data. What is shown: +complete GPS trace of the log file, positions of recorded signal +belonging to cell, little circles of different size indicating +signal strength (small circle: high signal, big circle: low signal), +and giant X the position of fitted tower location. + +--------------------------------------------------------- + +Operation: + +1) Set your logfile variable to the correct value +2) Set cellcheck = 1 (surely there is at least one cell in your log) +3) Set toplot = 1 +4) Run program - tells the total number of cells in logfile (i.e. the highest + cellcheck value one can use), and tries to fit the recorded signal + strengths. If there are not too many points or all at the same place, + don't bother, just discard that cell - until more data is logged. + Aim for the lowest RSS value (the "yyy" value in the RSS: xxx -> yyy line), + and run the program with the same cellcheck multiple times for that. +5) If you know that you have enough points for the cell you are looking at, + set toplot = 0 to get results quicker. +6) Collect fitted tower positions in datafile as: + mcc,mnc,cellid,lat,long +7) Ready to predict your position with "gsmpos" (hopefully) + + +...) If you know find better ways of predicting position, or bugs/errors in +the code, give feedback, please. Cheers! diff --git a/gsmpos/gsmfit.m b/gsmpos/gsmfit.m new file mode 100644 index 0000000..245dc6d --- /dev/null +++ b/gsmpos/gsmfit.m @@ -0,0 +1,204 @@ +%%%%%%%%%%%%%%%%%% +% (C) 2008 Gergely Imreh +% +% GPLv3 or later +%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%% +% Crude Octave3 code to attempt to find GSM tower positions +% using log from Cell Locator: +% http://repo.or.cz/w/CellLocator.git +% http://projects.openmoko.org/projects/rocinante/ +% +% Very rough, trying to fit free-space propagation model to +% recorded signal. Also visualizing log data. +%%%%%%%%%%%%%%%%%% + +function gsmfit() + +%%%%%%%%%%%% +% Parameters to modify +%%%%%%%%%%%% +% GSM/GPS Logfile +logfile = 'log_cells.csv'; +% Clean up before use! Remove headers, and possible invalid lines + +% Which GSM cell to check +cellcheck = 1; +% it is the number in the order of the cells listed in logfile +% cellcheck = [1..number-of-different-cellid] + +% Whether to plot: 0/1 +toplot = 1; + +%%%%End parameters + +disp('Parameters ==============') +disp(['Logfile: ',logfile]) +disp(['Cell to check : ' num2str(cellcheck)]) + +close all; +format none + +%% Load logfile +data = load(logfile); + + +%% Parsing log into separate variables +% GPS position +lon = data(:,3); +lat = data(:,4); +% CellID +cid = data(:,10); +% Base station id +bs = data(:,9); +% Received signal level +rxl = data(:,11); +% MCC/MNC numbers for the network +mcc = data(1,6); +mnc = data(1,7); + +%% Getting all distinct CellID numbers +cidlist = unique(cid); +% Sometimes there's unknown cellid logged as '0' - omit +if (cidlist(1) == 0) + cidlist = cidlist(2:end); +end % if +disp(['Number of cells in log: ' num2str(length(cidlist))]) +disp(['MCC,MNC: ' num2str(mcc) ',' num2str(mnc)]); + +if (cellcheck < 1) or (cellcheck > length(cidlist)) + disp('Invalid cell number - no such cell to check') +end % if + + +% plot complete logged GPS track +if toplot == 1 + hold off + figure(1) + plot(lon,lat,'g.') + axis equal + hold on +end % if + +out = []; + +%% Do fitting +%%%%% In loop to allow for later extending to batch procession +%%%%% Currently only one at a time: need of visual inspection of data for each cell +for k = cellcheck:cellcheck + + %% Current CellID + id = cidlist(k); + disp(['MCC,MNC,CellID: ' num2str(mcc) ',' num2str(mnc) ',' num2str(id)]) + + %% Number of logged points + nlog = length(cid(cid == id)); + disp(['Logged points: ' num2str(nlog)]) + + %% Select logged parameters for just this cell + clat = lat(cid == id); + clon = lon(cid == id); + crxl = rxl(cid == id); + cbs = bs(cid==id); + basestations = unique(cbs); + disp(["Base stations: " num2str(basestations')]) + + %% Turn RxL (dB) into signal scale + csig = 10.^(crxl/10); + + % Calculate all signal compared to maximum recorded + crxlnorm = max(crxl)-crxl; + % Distance scaling power: d^(-dsq) : dsq=2 is free space propagation + dsc = 2; + %%%% scale factor only for plotting, that should be handled more cleverly... + dscale = 0.00001; + + %% Approximate distance scale: + % RxL (dB) -> Signal (relative) -> Distance (with assumed scaling) + ds = 1./(10.^(crxlnorm/10)).^(1/dsc); + dist = dscale*1./ds; + + %% Plot data recorded for given cell + % Location helper: + % - circles around logged positions + % - circle radius is inversely proportional to signal strength + % ---> tower 'close': small circle; tower 'far': big circle + if toplot == 1 + figure(1) + plot(clon,clat,'b.') + for l = 1:nlog + c1 = dist(l)*cos(linspace(0,2*pi))+clon(l); + c2 = dist(l)*sin(linspace(0,2*pi))+clat(l); + plot(c1,c2) + end % for + axis equal + end % if + + + %% For fitting, guess the tower position close to the max recorded signal + %% an add random offset + [maxsig maxsignum] = max(crxl); + offsetm = 0.003; + fitofflon = offsetm*rand-offsetm/2; + fitofflat = offsetm*rand-offsetm/2; + b_in = [clon(maxsignum)+fitofflon clat(maxsignum)+fitofflat 0.01]; + + %% use the RxL for fitting... + [yhat, b_est] = leasqr([clon clat],crxl,b_in,"gsmsurf",0.00001,3000); + disp(['Fitted tower coordinates (Lat/Lon): ' num2str(b_est(2)) ',' num2str(b_est(1))]); + + % Residual sum of squares (for fit quality purposes: smaller the better! + RSSin = sum((gsmsurf([clon clat],b_in)-crxl).^2); + RSSest = sum((gsmsurf([clon clat],b_est)-crxl).^2); + disp(['RSS: ' num2str(RSSin) ' -> ' num2str(RSSest)]) + + %%% Put fitted tower position on map: big red X + if toplot == 1 + plot(b_est(1),b_est(2),'rx','MarkerSize',32); + xlabel(['CellID:' num2str(id) '/ RSS: ' num2str(RSSest) ' / Pos: ' num2str([b_est(1) b_est(2)])]) + axis equal; + hold off + end % if + + %% To assess fit quality: + if toplot == 1 + figure(2) + lon0 = b_est(1); + lat0 = b_est(2); + r = sqrt((clon-lon0).^2 + (clat-lat0).^2); +% yhat = yhat./log(10)*10; + plot(r,crxl,'@',r,yhat,'x') +% plot(r,lcsig,'@',r,yhat,'x') + ylabel('RxL (dB)') + xlabel(['Distance from fitted tower position / CellID:' num2str(id) ]) + end % if + +end % for +end % function gsmfit + + +%%%%%% +%%% Additional sub-routines +%%%%%% + +function yhat = gsmsurf(X,b) +%% Signal propagation surface according to propagation assumption + + %% Assumed distance scaling: d^(-dsc): dsc=2 is free-space propagation + dsc = 2; + %%%%% Separating input parameters into variables + % proposed cell position + lon0 = b(1); + lat0 = b(2); + % tower output scale factor - "signal at tower" - purely for fitting + s0 = b(3); + %%%%% + + % Distance from checked centre + r = sqrt((X(:,1)-lon0).^2 + (X(:,2)-lat0).^2); + % Estimated 'signal' + yhat = 10*log(s0*r.^(-dsc))/log(10); +end + +%%%%%%%%% END \ No newline at end of file -- 2.11.4.GIT