bluetooth le direction finding for tracking node position -凯发k8网页登录
this example shows you how to track the 2-d or 3-d position of a bluetooth® low energy (le) node by using bluetooth® toolbox.
using this example, you can:
simulate direction-finding packet exchange in the presence of radio frequency (rf) front end impairments, path loss, and additive white gaussian noise (awgn).
track the node position by using bluetooth direction-finding features and position estimation techniques.
improve the location accuracy by using a kalman filter from sensor fusion and tracking toolbox™.
bluetooth le direction finding
the bluetooth core specification v5.1 [2] introduced new features that support high-accuracy direction finding. bluetooth direction finding presents two distinct methods to estimate the 2-d or 3-d location of the bluetooth le node: angle of arrival (aoa) and angle of departure (aod). in aoa and aod techniques, the antenna array is present at the receiver and transmitter, respectively. some commonly used antenna arrays are uniform linear array (ula) and uniform rectangular array (ura). for more information on bluetooth direction finding, see .
location finding methods
the example uses these terminologies:
node - the device whose location is to be determined.
locator - the receiving device (in aoa calculation) or transmitting device (in aod calculation).
this figure shows two different ways to calculate the location of the bluetooth le node. each circle indicates the range of a bluetooth le locator.
distance plus angle: if one locator is present in the range of the bluetooth le node, use this method to estimate the node position. the location can be estimated by calculating the aoa or aod of a direction-finding signal and estimating the distance between the locator and node by performing a path loss calculation on the received signal strength indicator (rssi).
angulation: if two or more locators are present in the range of the bluetooth le node, use this method to estimate the node position. the location can be estimated by triangulating the aoa or aod between each locator and node.
simulation parameters
specify the linear motion model. generate dynamic node positions by using one of these motion models:
2-d constant velocity
3-d constant velocity
2-d constant acceleration
3-d constant acceleration
motionmodel ="2-d constant velocity"; % linear motion model
specify the direction-finding method, the direction-finding packet type, and the physical layer (phy) transmission mode.
dfmethod ="aoa"; % direction-finding method dfpackettype =
"connectioncte"; % direction-finding packet type % phy transmission mode, must be le1m or le2m (for connectioncte) and le1m % (for connectionlesscte) phymode =
"le1m";
specify the antenna array parameters.
% antenna array size, must be a scalar (represents ula) or a vector % (represents ura) for 2-d or 3-d positioning, respectively arraysize =16; % normalized spacing between the antenna elements with respect to % wavelength elementspacing =
0.5; % antenna switching pattern, must be a 1xm row vector and m must be in the % range [2, 74/slotduration 1] switchingpattern =
1:prod(arraysize);
specify the waveform generation and reception parameters.
slotduration =2; % slot duration in microseconds sps =
4; % samples per symbol channelindex =
17; % channel index crcinit =
'555551'; % cyclic redundancy check (crc) initialization accessaddress =
'01234567'; % access address payloadlength =
1; % payload length in bytes % length of constant tone extension (cte) in microseconds, must be in the % range [16, 160] with 8 microseconds step size ctelength =
160; % sample offset, must be a positive integer in the range [sps/8,7*sps/8] % for le1m and [sps/4,7*sps/4] for le2m sampleoffset =
2;
specify the bit energy to noise density ratio (ebno), environment and transmitter power.
ebno =25; % eb/no in db environment =
"outdoor"; % environment txpower =
0; % transmit power in dbm
validate the simulation parameters.
% configure the bluetoothrangeconfig object rangeconfig = bluetoothrangeconfig; rangeconfig.signalpowertype = "receivedsignalpower"; rangeconfig.environment = environment; rangeconfig.transmitterpower = txpower; rangeconfig.transmittercableloss = 0; % transmitter cable loss rangeconfig.receivercableloss = 0; % receiver cable loss numdimensions = 2 (strcmp(motionmodel,"3-d constant velocity") || strcmp(motionmodel,"3-d constant acceleration")); if numdimensions == 2 && size(arraysize,2) ~= 1 error("the arraysize must be a scalar for 2-d position estimation"); elseif numdimensions == 3 && size(arraysize,2) ~= 2 error("the arraysize must be a 1-by-2 vector for 3-d position estimation"); end if strcmp(dfpackettype,"connectioncte") && payloadlength ~= 1 error("the payloadlength must be 1 byte for direction-finding packet type of connectioncte"); elseif strcmp(dfpackettype,"connectionlesscte") && payloadlength < 3 error("the payloadlength must be greater than or equal to 3 bytes for direction-finding packet type of connectionlesscte"); end
create and configure a bluetooth le angle estimation configuration object.
if numdimensions == 3 && isscalar(elementspacing) elementspacing = [elementspacing elementspacing]; end cfg = bleangleestimateconfig("arraysize",arraysize, ... "slotduration",slotduration, ... "switchingpattern",switchingpattern, ... "elementspacing",elementspacing); validateconfig(cfg); % validate the configuration object pos = getelementposition(cfg); % element positions of an array % derive type of cte based on slot duration and direction-finding method if strcmp(dfmethod,"aoa") ctetype = [0;0]; else ctetype = [0;1]; if slotduration == 1 ctetype = [1;0]; end end
create and configure system objects™ for receiver processing.
% create and configure preamble detector system object accessaddbitslen = 32; accessaddbits = int2bit(hex2dec(accessaddress),accessaddbitslen,false); refsamples = helperblereferencewaveform(phymode,accessaddbits,sps); prbdet = comm.preambledetector(refsamples,"detections","first"); % create and configure coarse frequency compensator system object phyfactor = 1 strcmp(phymode,"le2m"); samplerate = 1e6*phyfactor*sps; coarsesync = comm.coarsefrequencycompensator("modulation","oqpsk",... "samplerate",samplerate,... "samplespersymbol",2*sps,... "frequencyresolution",100); % create and configure crc detector system object crclen = 24; % crc length crcdet = comm.crcdetector("x^24 x^10 x^9 x^6 x^4 x^3 x 1","directmethod",true,... "initialconditions",int2bit(hex2dec(crcinit),crclen).');
packet transmission and reception
perform these steps to track a bluetooth le node.
consider 15 locators and a node in a network. generate 30 node positions based on the motion of the node.
at each node position, consider active locators in 80 meters range.
model the direction-finding packet exchange between the node and each active locator. this figure shows the processing chain that the example uses to estimate the angle(s) and distances between the node and active locator.
4. repeat steps 2 and 3 for each node position.
5. estimate the node position by using the known active locator positions, distances, and angles.
rng('default'); % initialize the random number generator stream snr = ebno - 10*log10(sps); % signal to noise ratio (snr) in db headerlen = 16 8*strcmp(dfpackettype,"connectioncte"); % header length preamblelen = 8*phyfactor; % preamble length % derive initial state for whitening based on channel index dewhitenstatelen = 6; chanidxbin = int2bit(channelindex,dewhitenstatelen).'; initstate = [1 chanidxbin]; dewhiten = bluetoothwhiten(initialconditions=initstate'); % consider one node and 15 locators in a network. generate random locator % positions and generate node positions based on the motion model. numlocators = 15; % number of locators in a network numnodepositions = 30; % number of node positions in a network [posnode,poslocators] = helperblepositions(motionmodel,numnodepositions,numlocators); % simulate motion of the node by looping over the generated node positions. % for each node position, find active locators. posnodeest = zeros(numdimensions,numnodepositions); for inumnode = 1:numnodepositions % consider active locators that are within 80 meters of range to the % node [posactivelocators,angleactive,distanceactive] = helperbleactivelocators(posnode(:,inumnode),poslocators); poslocatorbuffer{:,inumnode} = posactivelocators; %#oknumactivelocators = size(posactivelocators,2); % number of active locators % estimate the pathloss between the transmitter and the receiver based % on the distance and the environment pllin = helperbluetoothestimatepathloss(rangeconfig.environment,distanceactive); angleest = zeros(numactivelocators,numdimensions-1); % estimated angles [pathlossdb,linkfailflag] = deal(zeros(numactivelocators,1)); % estimated path loss and flag to indicate link fail % loop over number of active locators for i = 1:numactivelocators
transmitter
% generate direction-finding packet data = helperblegeneratedfpdu(dfpackettype,ctelength,ctetype,payloadlength,crcinit); % generate bluetooth le waveform blewaveform = blewaveformgenerator(data,"mode",phymode,"samplespersymbol",sps,... "channelindex",channelindex,"dfpackettype",dfpackettype,"accessaddress",accessaddbits); % if the direction-finding method is aod, perform antenna switching if strcmp(dfmethod,"aod") blewaveform = helperblesteerswitchantenna(blewaveform,angleactive(i,:),... phymode,sps,dfpackettype,payloadlength,cfg); end
channel
add rf impairments to the waveform. attenuate the impaired waveform according to the transmitter power and pathloss.
% add frequency offset freqoffset = randsrc(1,1,-10e3:100:10e3); freqwaveform = helperblefrequencyoffset(blewaveform,samplerate,freqoffset); % add timing offset timingoff = randsrc(1,1,1:0.2:20); timingoffwaveform = helperbledelaysignal(freqwaveform,timingoff); % add dc offset dcvalue = (5/100)*max(timingoffwaveform); dcwaveform = timingoffwaveform dcvalue; % initialize phase noise system object phasenoise = comm.phasenoise(level=[-130 -136],frequencyoffset=[1e4 1e5],samplerate=samplerate); % add phase noise noisywaveform = phasenoise(dcwaveform); % release the phase noise object release(phasenoise); % apply transmitter power and path loss dbmconverter = 30; tximpairedwaveform = 10^((rangeconfig.transmitterpower-dbmconverter)/20)*noisywaveform/pllin(i);
receiver
perform these steps to recover iq samples for angle estimation.
perform antenna switching (for aoa)
perform packet detection
recover the header bits
perform crc detection
perform cte iq sampling
% perform antenna steering if the direction-finding method is aoa if strcmp(dfmethod,"aoa") steervec = helperblesteeringvector(angleactive(i,:),pos); steeredwaveform = tximpairedwaveform .* steervec.'; else steeredwaveform = tximpairedwaveform; end % initialize the variables used for receiver processing packetdetect = 0; % packet detect flag headerflag = 1; % header flag to perform header decoding pducrcflag = 0; % pdu and crc flag to perform crc detection samplingflag = 0; % cte iq sampling flag crcerror = 1; % checksum flag samplesperframe = 8*sps; % samples considered for packet detection samplespermodule = samplesperframe; % samples considered for header, pdu and crc, cte iq sampling numtimes = ceil(length(timingoffwaveform)/samplesperframe) 1; % divide the waveform as per the samples considered countpacketdetect = 0; % counter for packet detection modulestartind = 0; % start index for each module % append zeros to the waveform rxwaveformbuffer = [steeredwaveform;zeros(samplesperframe*numtimes-length(steeredwaveform) samplesperframe,size(steeredwaveform,2))]; for j = 1:numtimes % retrieve the waveform corresponding to each module if (j-1)*samplesperframe modulestartind samplespermodule <= length(rxwaveformbuffer) rxchunkwaveform = rxwaveformbuffer((j-1)*samplesperframe modulestartind (1:samplespermodule),:); else rxchunkwaveform = rxwaveformbuffer((j-1)*samplesperframe modulestartind 1:end,:); end
a sampling flag is used to handle antenna switching for aoa.
if the sampling flag is disabled, use the waveform from the first antenna. add awgn to the waveform.
if the sampling flag is enabled, perform antenna switching (for aoa). add awgn to the waveform. perform iq sampling on the cte samples.
if ~samplingflag rxnoisywaveform = awgn(rxchunkwaveform(:,1),snr,"measured"); % add awgn else if strcmp(dfmethod,"aoa") % perform antenna switching rxswitchwaveform = helperbleswitchantenna(rxchunkwaveform,phymode,sps,slotduration,switchingpattern); else rxswitchwaveform = rxchunkwaveform; end rxnoisywaveform = awgn(rxswitchwaveform,snr,"measured"); % add awgn % perform iq sampling on the cte ctesamples = rxnoisywaveform(1:ctetime*8*sps*phyfactor); iqsamples = blecteiqsample(ctesamples,"mode",phymode,... "samplespersymbol",sps,"slotduration",slotduration,"sampleoffset",sampleoffset); samplingflag = 0; break; end
perform these steps for packet detection:
load the waveform into a buffer to perform preamble detection.
estimate and compensate dc and frequency offsets.
perform gaussian minimum shift keying (gmsk) demodulation.
compare the decoded access address with the known access address. if these addresses match, the packet is detected.
if packetdetect == 0 countpacketdetect = countpacketdetect 1; rcvsigbuffer((countpacketdetect-1)*samplesperframe 1:(countpacketdetect)*samplesperframe) = rxnoisywaveform; if countpacketdetect >= (preamblelen accessaddbitslen headerlen)*sps/samplesperframe % perform timing synchronization [~, dtmt] = prbdet(rcvsigbuffer.'); release(prbdet) prbdet.threshold = max(dtmt); prbidx = prbdet(rcvsigbuffer.'); release(prbdet) if prbidx >=length(refsamples) rcvtrim = rcvsigbuffer(1 prbidx-length(refsamples):end).'; else rcvtrim = rcvsigbuffer.'; end % estimate and compensate dc offset estimateddcoffset = mean(rcvtrim(1:length(refsamples)))-mean(refsamples)*sqrt(var(rcvtrim)); rcvdcfree = rcvtrim-estimateddcoffset; % estimate and compensate frequency offset [~,freqoff] = coarsesync(rcvdcfree); release(coarsesync) [rcvfreqoffsetfree,inifreqstate] = helperblefrequencyoffset(rcvdcfree,samplerate,-freqoff); % perform gmsk demodulation x = rem(length(rcvfreqoffsetfree),sps); if x remsad = sps - x; else remsad = x; end remnumsamples = x; remsamples = rcvfreqoffsetfree(end-remnumsamples 1:end); [demodsoftbits,demodinitphase] = helperblegmskdemod(rcvfreqoffsetfree(1:end-remnumsamples),sps,0); demodbits = demodsoftbits>0; % check for access address match if length(demodbits) >= preamblelen accessaddbitslen accessaddress = int8(demodbits(preamblelen (1:accessaddbitslen))>0); decodedata = int8(demodbits(accessaddbitslen preamblelen 1:end)>0); if isequal(accessaddbits,accessaddress) packetdetect = 1; samplespermodule = headerlen*sps-length(decodedata)*sps remsad; rcvsigbuffer = []; end end end end
recover the header bits if the packet is detected.
if packetdetect if headerflag && (length(rxnoisywaveform) ~= samplesperframe || samplespermodule <= 0) % remove dc offset rxdcfreewaveform = rxnoisywaveform-estimateddcoffset; % frequency offset correction [rxnoisywaveform,inifreqstate] = helperblefrequencyoffset(rxdcfreewaveform,samplerate,-freqoff,inifreqstate); % perform gmsk demodulation on the header bits if (length(rxnoisywaveform) ~= samplesperframe) rcvsigheaderbuffer = [remsamples;rxnoisywaveform]; [headersoftbits,demodinitphase] = helperblegmskdemod(rcvsigheaderbuffer,sps,demodinitphase); decodedataheader = [decodedata;headersoftbits>0]; elseif samplespermodule <= 0 decodedataheader = decodedata; end % perform data dewhitening dewhitenedbits = dewhiten(decodedataheader); % extract pdu length pdulenfield = double(dewhitenedbits(9:16)); % second byte of pdu header pduleninbits = bi2de(pdulenfield')*8; % retrieve cte time and cte type information from bits if strcmp(dfpackettype,"connectioncte") [ctetime,ctetypedec] = helperblecteinfoextract(dewhitenedbits,dfpackettype); % slot duration based on cte type if ctetypedec == 1 || ctetypedec == 2 slotduration = ctetypedec; else slotduration = cfg.slotduration; end end % update the samples per frame, start index to perform % crc detection if samplespermodule modulestartind = samplespermodule-samplesperframe; else modulestartind = 0; end samplespermodule = (pduleninbits crclen-length(dewhitenedbits) headerlen)*sps; pducrcflag = 1; % enable the flag for crc detection headerflag = 0; % disable the header recovery flag rcvsigheadervar = rxnoisywaveform; % load the waveform for pathloss estimation rxnoisywaveform = []; end
perform crc detection based on the pdu length recovered from header bits.
if pducrcflag && ~isempty(rxnoisywaveform) % remove dc offset rxdcfreewaveform = rxnoisywaveform-estimateddcoffset; % frequency offset correction [rxnoisywaveform,inifreqstate] = helperblefrequencyoffset(rxdcfreewaveform,samplerate,-freqoff,inifreqstate); % perform gmsk demodulation on the pdu and crc bits x = rem(length(rxnoisywaveform),sps); if x remnumsamples = sps - x; else remnumsamples = x; end demodpducrc = helperblegmskdemod([rxnoisywaveform;zeros(remnumsamples,1)],sps,demodinitphase)>0; % perform data dewhitening dewhitenedpducrc = dewhiten(demodpducrc); reset(dewhiten) % retrieve cte time and cte type information from bits headerpducrc = [double(dewhitenedbits);dewhitenedpducrc]; if strcmp(dfpackettype,"connectionlesscte") [ctetime,ctetypedec] = helperblecteinfoextract(headerpducrc,dfpackettype); % slot duration based on cte type if ctetypedec == 1 || ctetypedec == 2 slotduration = ctetypedec; else slotduration = cfg.slotduration; end end % perform crc detection [dfpdu,crcerror] = crcdet(headerpducrc); if crcerror break; end % update the samples per frame, start index to perform % iq sampling modulestartind = samplespermodule-samplesperframe modulestartind; samplespermodule = ctetime*8*sps*phyfactor; samplingflag = 1; % enable the flag to perform iq sampling pducrcflag = 0; % disable the pdu and crc flag rcvsigpduvar = rxnoisywaveform; % load the waveform for pathloss estimation end end end
estimate angle and pathloss
use the iq samples from the receiver to estimate the angles and path loss.
perform angle estimation if both of these conditions meet:
crc is detected.
a minimum number of non-zero iq samples are present in
iqsamples.
% iq samples must contain at least eight samples from reference % period and one sample from each antenna refsamplelength = 8; % reference samples length miniqsamples = refsamplelength getnumelements(cfg)-1; if ~crcerror && length(nonzeros(iqsamples)) >= miniqsamples % if packet detection is successful angleest(i,:) = bleangleestimate(iqsamples,cfg); rangeconfig.receivedsignalpower = 10*log10(var([rcvfreqoffsetfree;rcvsigheadervar;rcvsigpduvar])) 30; pathlossdb(i) = pathloss(rangeconfig); else linkfailflag(i) = 1; % if packet detection fails, enable the link fail flag end end
estimate node position
estimate the position of the node in the network by using one of these two methods.
angulation: if two or more number of active locators are successfully decoded.
distance-angle: if one active locator is successfully decoded.
if any(linkfailflag == 1) idx = find(linkfailflag == 1); posactivelocators(:,idx) = []; angleest(idx,:) = []; pathlossdb(idx) = []; end if (numactivelocators-nnz(linkfailflag)) >= 2 && ~all(angleest(:,1) == angleest(1,1)) % estimate the node position by using the angulation method posnodeest(:,inumnode) = blepositionestimate(posactivelocators,"angulation",angleest.'); elseif (numactivelocators-nnz(linkfailflag)) == 1 || (~all(linkfailflag == 1) && all(angleest(:,1) == angleest(1,1))) % estimate the distance between the transmitter and the receiver % based on the environment range = bluetoothrange(rangeconfig); distancevalues = min(range):max(range); idx = randi(length(distancevalues),1); distanceest = distancevalues(idx) % estimate the node position by using the distance-angle method posnodeest(:,inumnode) = blepositionestimate(posactivelocators,"distance-angle",distanceest.',angleest.'); else posnodeest(:,inumnode) = nan(numdimensions,1); end end
simulation results
visualize the bluetooth le node tracking.
poserr = sqrt(sum((posnodeest-posnode).^2)); % compute the position error disp(["positioning error in meters = ", num2str(poserr)])
"positioning error in meters = " "0.1423 0.3255 0.14592 0.525 0.25757 0.64587 0.11686 0.27105 0.19259 0.11363 0.20299 0.20931 0.10998 0.53644 0.84425 0.70373 0.2665 0.17448 0.95019 0.37859 0.12991 0.7048 0.5538 0.11129 0.064095 0.68762 0.58582 0.083517 0.045576 0.089482"
if ~all(isnan(posnodeest(1,:))) helperblevisualizenodetracking(poslocators,posnode,poslocatorbuffer,posnodeest) end
this example enables you to estimate the 2-d or 3-d position of a bluetooth le node using the direction-finding functionality and location estimation techniques. the example also measures positioning accuracy at each node position.
further exploration
improve the position accuracy of the bluetooth le node by using the kalman filter from sensor fusion and tracking toolbox™. the kalman filter is a recursive algorithm that estimates the track of an object by updating the current state using the previous state. when the evolution of the state follows a linear motion model and the measurements are linear functions of the state, the kalman filter is linear. for more information about the linear kalman filter, see (sensor fusion and tracking toolbox).
the two important characteristics of kalman filter are:
predicting the node's future location
reducing the noise introduced by inaccurate detections
if node estimation measurement is available, the kalman filter reduces the noise by predicting and correcting together. if node estimation measurement is not available (all the node-locators links fail), the kalman filter estimates the location based on the previous measurements. as there is no prior knowledge on initial estimates of the velocity or acceleration of node, the position estimation can be transient during the initial phase. use this sample code snippet after estimating node positions.
% if all(isnan(posnodeest(1,:))) % disp("all direction-finding packet transmissions failed") % else % posnodetrackest = helperblekalmanfilter(motionmodel,posnodeest); % postrackerr = sqrt(sum((posnodetrackest-posnode).^2)); % compute the position error with kalman filter % disp(["positioning error with kalman filter in meters = ",num2str(postrackerr)]) % end
appendix
the example uses these helpers:
: generate bluetooth le reference samples
: generate locators and node positions
: generate active locators parameters
: estimate the path loss between the node and locator
: generate direction-finding packet pdu
: perform antenna steering and switching
: apply frequency offset to the input signal
: introduce time delay in the signal
: generate steering vector
: perform antenna switching
: perform gmsk demodulation
: extract the cte information
: track the bluetooth le node using linear kalman filter
: generate bluetooth le node tracking visualization
selected bibliography
“bluetooth® technology website – the official website for the bluetooth wireless technology. get up to date specifications, news, and development info.” accessed january 4, 2023.
“core specification – bluetooth® technology website.” accessed january 4, 2023.
“core specification – bluetooth® technology website.” accessed january 4, 2023. .
see also
functions
- | | | |
objects
- |
related topics
select a web site
choose a web site to get translated content where available and see local events and offers. based on your location, we recommend that you select: .
you can also select a web site from the following list:
americas
- (español)
- (english)
- (english)
europe
- (english)
- (english)
- (deutsch)
- (español)
- (english)
- (français)
- (english)
- (italiano)
- (english)
- (english)
- (english)
- (deutsch)
- (english)
- (english)
- switzerland
- (english)
asia pacific
- (english)
- (english)
- (english)
- 中国
- (日本語)
- (한국어)