%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % \file LoRaPHY.m % % \brief Physical Layer LoRa Modulator/Demodulator/Encoder/Decoder % % \version 0.2.1 % % \repo https://github.com/jkadbear/LoRaPHY % % \copyright MIT License, 2020-2022 % % \author jkadbear %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% classdef LoRaPHY < handle & matlab.mixin.Copyable %LORAPHY LoRa physical layer implementation %%% Example %%% % rf_freq = 470e6; % sf = 12; % bw = 125e3; % fs = 1e6; % phy = LoRaPHY(rf_freq, sf, bw, fs); % phy.has_header = 1; % explicit header mode % symbols = [2541,1153,673,2397,1189,3509,41,3089,3237,3917,2729,2765,1417,2833,1389,801,3197,345,961,745,3101,297,1893,469]'; % [data, checksum] = phy.decode(symbols); % disp(data); % CODE: 09 90 40 01 02 03 04 05 06 07 08 09 BA 2E % disp(checksum); properties rf_freq % carrier frequency sf % spreading factor (7,8,9,10,11,12) bw % bandwidth (125kHz 250kHz 500kHz) fs % sampling frequency cr % code rate: (1:4/5 2:4/6 3:4/7 4:4/8) payload_len % payload length has_header % explicit header: 1, implicit header: 0 crc % crc = 1 if CRC Check is enabled else 0 ldr % ldr = 1 if Low Data Rate Optimization is enabled else 0 whitening_seq % whitening sequence crc_generator % CRC generator with polynomial x^16+x^12+x^5+1 header_checksum_matrix % we use a 12 x 5 matrix to calculate header checksum preamble_len % preamble length sig % input baseband signal downchirp % ideal chirp with decreasing frequency from B/2 to -B/2 upchirp % ideal chirp with increasing frequency from -B/2 to B/2 sample_num % number of sample points per symbol bin_num % number of bins after FFT (with zero padding) zero_padding_ratio % FFT zero padding ratio fft_len % FFT size preamble_bin % reference bin in current decoding window, used to eliminate CFO cfo % carrier frequency offset fast_mode % set `true` for fast execution (ignore low-pass filter) is_debug % set `true` for debug information hamming_decoding_en % enable hamming decoding end methods function self = LoRaPHY(rf_freq, sf, bw, fs) %LORAPHY Construct an instance of this class % Hexadecimal or binary values representation require at least % MATLAB R2019b % https://www.mathworks.com/help/matlab/matlab_prog/specify-hexadecimal-and-binary-numbers.html if verLessThan('matlab', '9.7') error('Error. Newer version of MATLAB is required ( >=R2019b ).'); end self.rf_freq = rf_freq; self.sf = sf; self.bw = bw; self.fs = fs; self.has_header = 1; self.crc = 1; self.fast_mode = false; self.is_debug = false; self.hamming_decoding_en = true; self.zero_padding_ratio = 10; self.cfo = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The whitening sequence is generated by an LFSR % x^8+x^6+x^5+x^4+1 % Use the code below to generate such sequence %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % reg = 0xFF; % for i = 1:255 % fprintf("0x%x, ", reg); % reg = bitxor(bitshift(reg,1), bitxor(bitget(reg,8), bitxor(bitget(reg,6), bitxor(bitget(reg,5), bitget(reg,4))))); % end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% self.whitening_seq = uint8([0xff, 0xfe, 0xfc, 0xf8, 0xf0, 0xe1, 0xc2, 0x85, 0xb, 0x17, 0x2f, 0x5e, 0xbc, 0x78, 0xf1, 0xe3, 0xc6, 0x8d, 0x1a, 0x34, 0x68, 0xd0, 0xa0, 0x40, 0x80, 0x1, 0x2, 0x4, 0x8, 0x11, 0x23, 0x47, 0x8e, 0x1c, 0x38, 0x71, 0xe2, 0xc4, 0x89, 0x12, 0x25, 0x4b, 0x97, 0x2e, 0x5c, 0xb8, 0x70, 0xe0, 0xc0, 0x81, 0x3, 0x6, 0xc, 0x19, 0x32, 0x64, 0xc9, 0x92, 0x24, 0x49, 0x93, 0x26, 0x4d, 0x9b, 0x37, 0x6e, 0xdc, 0xb9, 0x72, 0xe4, 0xc8, 0x90, 0x20, 0x41, 0x82, 0x5, 0xa, 0x15, 0x2b, 0x56, 0xad, 0x5b, 0xb6, 0x6d, 0xda, 0xb5, 0x6b, 0xd6, 0xac, 0x59, 0xb2, 0x65, 0xcb, 0x96, 0x2c, 0x58, 0xb0, 0x61, 0xc3, 0x87, 0xf, 0x1f, 0x3e, 0x7d, 0xfb, 0xf6, 0xed, 0xdb, 0xb7, 0x6f, 0xde, 0xbd, 0x7a, 0xf5, 0xeb, 0xd7, 0xae, 0x5d, 0xba, 0x74, 0xe8, 0xd1, 0xa2, 0x44, 0x88, 0x10, 0x21, 0x43, 0x86, 0xd, 0x1b, 0x36, 0x6c, 0xd8, 0xb1, 0x63, 0xc7, 0x8f, 0x1e, 0x3c, 0x79, 0xf3, 0xe7, 0xce, 0x9c, 0x39, 0x73, 0xe6, 0xcc, 0x98, 0x31, 0x62, 0xc5, 0x8b, 0x16, 0x2d, 0x5a, 0xb4, 0x69, 0xd2, 0xa4, 0x48, 0x91, 0x22, 0x45, 0x8a, 0x14, 0x29, 0x52, 0xa5, 0x4a, 0x95, 0x2a, 0x54, 0xa9, 0x53, 0xa7, 0x4e, 0x9d, 0x3b, 0x77, 0xee, 0xdd, 0xbb, 0x76, 0xec, 0xd9, 0xb3, 0x67, 0xcf, 0x9e, 0x3d, 0x7b, 0xf7, 0xef, 0xdf, 0xbf, 0x7e, 0xfd, 0xfa, 0xf4, 0xe9, 0xd3, 0xa6, 0x4c, 0x99, 0x33, 0x66, 0xcd, 0x9a, 0x35, 0x6a, 0xd4, 0xa8, 0x51, 0xa3, 0x46, 0x8c, 0x18, 0x30, 0x60, 0xc1, 0x83, 0x7, 0xe, 0x1d, 0x3a, 0x75, 0xea, 0xd5, 0xaa, 0x55, 0xab, 0x57, 0xaf, 0x5f, 0xbe, 0x7c, 0xf9, 0xf2, 0xe5, 0xca, 0x94, 0x28, 0x50, 0xa1, 0x42, 0x84, 0x9, 0x13, 0x27, 0x4f, 0x9f, 0x3f, 0x7f]'); self.header_checksum_matrix = gf([ 1 1 1 1 0 0 0 0 0 0 0 0 1 0 0 0 1 1 1 0 0 0 0 1 0 1 0 0 1 0 0 1 1 0 1 0 0 0 1 0 0 1 0 1 0 1 1 1 0 0 0 1 0 0 1 0 1 1 1 1 ]); self.crc_generator = comm.CRCGenerator('Polynomial','X^16 + X^12 + X^5 + 1'); self.preamble_len = 6; self.init(); end function init(self) % init Initialize some parameters self.bin_num = 2^self.sf*self.zero_padding_ratio; self.sample_num = 2*2^self.sf; self.fft_len = self.sample_num*self.zero_padding_ratio; self.downchirp = LoRaPHY.chirp(false, self.sf, self.bw, 2*self.bw, 0, self.cfo, 0); self.upchirp = LoRaPHY.chirp(true, self.sf, self.bw, 2*self.bw, 0, self.cfo, 0); % Low Data Rate Optimization (LDRO) mode in LoRa % If the chirp peird is larger than 16ms, the least significant % two bits are considered unreliable and are neglected. if 2^(self.sf)/self.bw > 16e-3 self.ldr = 1; else self.ldr = 0; end end function pk = dechirp(self, x, is_up) % dechirp Apply dechirping on the symbol starts from index x % % input: % x: Start index of a symbol % is_up: `true` if applying up-chirp dechirping % `false` if applying down-chirp dechirping % output: % pk: Peak in FFT results of dechirping % pk = (height, index) if nargin == 3 && ~is_up c = self.upchirp; else c = self.downchirp; end ft = fft(self.sig(x:x+self.sample_num-1).*c, self.fft_len); ft_ = abs(ft(1:self.bin_num)) + abs(ft(self.fft_len-self.bin_num+1:self.fft_len)); pk = LoRaPHY.topn([ft_ (1:self.bin_num).'], 1); end function x = detect(self, start_idx) % detect Detect preamble % % input: % start_idx: Start index for detection % output: % x: Before index x, a preamble is detected. % x = -1 if no preamble detected ii = start_idx; pk_bin_list = []; % preamble peak bin list while ii < length(self.sig)-self.sample_num*self.preamble_len % search preamble_len-1 basic upchirps if length(pk_bin_list) == self.preamble_len - 1 % preamble detected % coarse alignment: first shift the up peak to position 0 % current sampling frequency = 2 * bandwidth x = ii - round((pk_bin_list(end)-1)/self.zero_padding_ratio*2); return; end pk0 = self.dechirp(ii); if ~isempty(pk_bin_list) bin_diff = mod(pk_bin_list(end)-pk0(2), self.bin_num); if bin_diff > self.bin_num/2 bin_diff = self.bin_num - bin_diff; end if bin_diff <= self.zero_padding_ratio pk_bin_list = [pk_bin_list; pk0(2)]; else pk_bin_list = pk0(2); end else pk_bin_list = pk0(2); end ii = ii + self.sample_num; end x = -1; end function [symbols_m, cfo_m, netid_m] = demodulate(self, sig) % demodulate LoRa packet demodulation % % input: % sig: Baseband signal in complex % output: % symbols_m: A matrix containing the demodulated results. % Each column vector represents the symbols of % a successfully demodulated packet. % cfo_m: A vector containing the carrier frequency offset % results. Each element represents the CFO of the % packet in symbols_m. self.cfo = 0; self.init(); if ~self.fast_mode sig = lowpass(sig, self.bw/2, self.fs); end % resample signal with 2*bandwidth self.sig = resample(sig, 2*self.bw, self.fs); symbols_m = []; cfo_m = []; netid_m = []; x = 1; while x < length(self.sig) x = self.detect(x); if x < 0 break; end % align symbols with SFD x = self.sync(x); % NetID pk_netid1 = self.dechirp(round(x-4.25*self.sample_num)); pk_netid2 = self.dechirp(round(x-3.25*self.sample_num)); netid_m = [netid_m; [mod((pk_netid1(2)+self.bin_num-self.preamble_bin)/self.zero_padding_ratio, 2^self.sf), ... mod((pk_netid2(2)+self.bin_num-self.preamble_bin)/self.zero_padding_ratio, 2^self.sf)] ]; % the goal is to extract payload_len from PHY header % header is in the first 8 symbols symbols = []; pk_list = []; if x > length(self.sig) - 8*self.sample_num + 1 return; end for ii = 0:7 pk = self.dechirp(x+ii*self.sample_num); pk_list = [pk_list; pk]; symbols = [symbols; mod((pk(2)+self.bin_num-self.preamble_bin)/self.zero_padding_ratio, 2^self.sf)]; end if self.has_header is_valid = self.parse_header(symbols); if ~is_valid x = x + 7*self.sample_num; continue; end end % number of symbols in the packet sym_num = self.calc_sym_num(self.payload_len); % demodulate the rest LoRa data symbols if x > length(self.sig) - sym_num*self.sample_num + 1 return; end for ii = 8:sym_num-1 pk = self.dechirp(x+ii*self.sample_num); pk_list = [pk_list; pk]; symbols = [symbols; mod((pk(2)+self.bin_num-self.preamble_bin)/self.zero_padding_ratio, 2^self.sf)]; end x = x + sym_num*self.sample_num; % compensate CFO drift symbols = self.dynamic_compensation(symbols); symbols_m = [symbols_m mod(round(symbols),2^self.sf)]; cfo_m = [cfo_m self.cfo]; end if isempty(symbols_m) warning('No preamble detected!'); end end function is_valid = parse_header(self, data) % parse_header Parse LoRa PHY header and set parameters % % input: % data: An eight elements vector containing header symbols % output: % is_valid: `true` if the header is valid % `false` if the header is invalid % compensate CFO drift symbols = self.dynamic_compensation(data); % gray coding symbols_g = self.gray_coding(symbols); % deinterleave codewords = self.diag_deinterleave(symbols_g(1:8), self.sf-2); % parse header nibbles = self.hamming_decode(codewords, 8); self.payload_len = double(nibbles(1)*16 + nibbles(2)); self.crc = double(bitand(nibbles(3), 1)); self.cr = double(bitshift(nibbles(3), -1)); % we only calculate header checksum on the first three nibbles % the valid header checksum is considered to be 5 bits % other 3 bits require further reverse engineering header_checksum = [bitand(nibbles(4), 1); de2bi(nibbles(5), 4, 'left-msb')']; header_checksum_calc = self.header_checksum_matrix * gf(reshape(de2bi(nibbles(1:3), 4, 'left-msb')', [], 1)); if any(header_checksum ~= header_checksum_calc) warning('Invalid header checksum!'); is_valid = 0; else is_valid = 1; end end function s = modulate(self, symbols) % modulate Modulate a baseband signal % % input: % symbols: A vector of chirp symbols to be modulated % valid symbol range: 0 to 2^sf-1 % output: % s: A valid LoRa baseband signal uc = LoRaPHY.chirp(true, self.sf, self.bw, self.fs, 0, self.cfo, 0); dc = LoRaPHY.chirp(false, self.sf, self.bw, self.fs, 0, self.cfo, 0); preamble = repmat(uc, self.preamble_len, 1); netid = [LoRaPHY.chirp(true, self.sf, self.bw, self.fs, 24, self.cfo, 0); LoRaPHY.chirp(true, self.sf, self.bw, self.fs, 32, self.cfo, 0)]; chirp_len = length(uc); sfd = [dc; dc; dc(1:round(chirp_len/4))]; data = zeros(length(symbols)*chirp_len, 1); for i = 1:length(symbols) data((i-1)*chirp_len+1:i*chirp_len) = LoRaPHY.chirp(true, self.sf, self.bw, self.fs, symbols(i), self.cfo, 0); end s = [preamble; netid; sfd; data]; end function symbols = encode(self, payload) % encode Encode bytes to symbols % % input: % payload: Payload of LoRa packet % output: % symbols: A vector representing the symbols of the packet if self.crc data = uint8([payload; self.calc_crc(payload)]); else data = uint8(payload); end plen = length(payload); sym_num = self.calc_sym_num(plen); % filling all symbols needs nibble_num nibbles nibble_num = self.sf - 2 + (sym_num-8)/(self.cr+4)*(self.sf-2*self.ldr); data_w = uint8([data; 255*ones(ceil((nibble_num-2*length(data))/2), 1)]); data_w(1:plen) = self.whiten(data_w(1:plen)); data_nibbles = uint8(zeros(nibble_num, 1)); for i = 1:nibble_num idx = ceil(i/2); if mod(i, 2) == 1 data_nibbles(i) = bitand(data_w(idx), 0xf); else data_nibbles(i) = bitshift(data_w(idx), -4); end end if self.has_header header_nibbles = self.gen_header(plen); else header_nibbles = []; end codewords = self.hamming_encode([header_nibbles; data_nibbles]); % interleave % first 8 symbols use CR=4/8 symbols_i = self.diag_interleave(codewords(1:self.sf-2), 8); ppm = self.sf - 2*self.ldr; rdd = self.cr + 4; for i = self.sf-1:ppm:length(codewords)-ppm+1 symbols_i = [symbols_i; self.diag_interleave(codewords(i:i+ppm-1), rdd)]; end symbols = self.gray_decoding(symbols_i); end function header_nibbles = gen_header(self, plen) % gen_header Generate a valid LoRa header % % input: % plen: Payload length % output: % header_nibbles: Header in nibbles header_nibbles = zeros(5, 1); header_nibbles(1) = bitshift(plen, -4); header_nibbles(2) = bitand(plen, 15); header_nibbles(3) = bitor(2*self.cr, self.crc); header_checksum = self.header_checksum_matrix * gf(reshape(de2bi(header_nibbles(1:3), 4, 'left-msb')', [], 1)); x = header_checksum.x; header_nibbles(4) = x(1); for i = 1:4 header_nibbles(5) = bitor(header_nibbles(5), x(i+1)*2^(4-i)); end end function checksum = calc_crc(self, data) % calc_crc Calculate payload CRC % % input: % data: Data in bytes % output: % checksum: CRC result switch length(data) case 0 checksum = [0; 0]; case 1 checksum = [data(end); 0]; case 2 checksum = [data(end); data(end-1)]; otherwise input = data(1:end-2); seq = self.crc_generator(reshape(logical(de2bi(input, 8, 'left-msb'))', [], 1)); checksum_b1 = bitxor(bi2de(seq(end-7:end)', 'left-msb'), data(end)); checksum_b2 = bitxor(bi2de(seq(end-15:end-8)', 'left-msb'), data(end-1)); checksum = [checksum_b1; checksum_b2]; end end function data_w = whiten(self, data) % whiten Whitening process in LoRa % % input: % data: Data in bytes % output: % data_w: Data after whitening len = length(data); data_w = bitxor(data(1:len), self.whitening_seq(1:len)); self.print_bin("Whiten", data_w); end function codewords = hamming_encode(self, nibbles) % hamming_encode Hamming encoding process in LoRa % % input: % nibbles: Data in nibbles % output: % codewords: Data after hamming encoding nibble_num = length(nibbles); codewords = uint8(zeros(nibble_num, 1)); for i = 1:nibble_num nibble = nibbles(i); p1 = LoRaPHY.bit_reduce(@bitxor, nibble, [1 3 4]); p2 = LoRaPHY.bit_reduce(@bitxor, nibble, [1 2 4]); p3 = LoRaPHY.bit_reduce(@bitxor, nibble, [1 2 3]); p4 = LoRaPHY.bit_reduce(@bitxor, nibble, [1 2 3 4]); p5 = LoRaPHY.bit_reduce(@bitxor, nibble, [2 3 4]); if i <= self.sf - 2 % the first SF-2 nibbles use CR=4/8 cr_now = 4; else cr_now = self.cr; end switch cr_now case 1 codewords(i) = bitor(uint8(p4)*16, nibble); case 2 codewords(i) = LoRaPHY.word_reduce(@bitor, [uint8(p5)*32 uint8(p3)*16 nibble]); case 3 codewords(i) = LoRaPHY.word_reduce(@bitor, [uint8(p2)*64 uint8(p5)*32 uint8(p3)*16 nibble]); case 4 codewords(i) = LoRaPHY.word_reduce(@bitor, [uint8(p1)*128 uint8(p2)*64 uint8(p5)*32 uint8(p3)*16 nibble]); otherwise % THIS CASE SHOULD NOT HAPPEN error('Invalid Code Rate!'); end end end function symbols_i = diag_interleave(self, codewords, rdd) % diag_interleave Diagonal interleaving % % input: % codewords: Data in nibbles % rdd: Bits with redundancy % For example, code rate 4/5 means rdd = 5 % output: % symbols_i: Symbols after diagonal interleaving tmp = de2bi(codewords, rdd, 'right-msb'); symbols_i = uint16(bi2de(cell2mat(arrayfun(@(x) circshift(tmp(:,x), 1-x), 1:rdd, 'un', 0))')); self.print_bin("Interleave", symbols_i); end function symbols = gray_decoding(self, symbols_i) % gray_decoding Gray decoding % `gray_decoding` is used in the ENCODING process % % input: % symbols_i: Interleaved symbols % output: % symbols: Final symbols to be modulated in a packet symbols = zeros(length(symbols_i), 1); for i = 1:length(symbols_i) num = uint16(symbols_i(i)); mask = bitshift(num, -1); while mask ~= 0 num = bitxor(num, mask); mask = bitshift(mask, -1); end if i <= 8 || self.ldr symbols(i) = mod(num * 4 + 1, 2^self.sf); else symbols(i) = mod(num + 1, 2^self.sf); end end end function sym_num = calc_sym_num(self, plen) % calc_sym_num Calculate number of symbols % % input: % plen: Payload length % output: % sym_num: Number of symbols sym_num = double(8 + max((4+self.cr)*ceil(double((2*plen-self.sf+7+4*self.crc-5*(1-self.has_header)))/double(self.sf-2*self.ldr)), 0)); end function plen = calc_payload_len(self, slen, no_redundant_bytes) % calc_payload_len Calculate payload length % % input: % slen: Number of symbols % no_redundant_bytes: `true` to fill the 0.5 redundant bytes % output: % plen: Payload length if nargin < 3 no_redundant_bytes = false; end % plen_float possibly has fractional part 0.5, which means % there would be 0.5 uncontrollable redundant byte in a packet. % The 0.5 byte results in unexpected symbols when called by % function `symbols_to_bytes`. To make all specified symbols % controllable, we use `ceil` instead of `floor` when % no_redundant_bytes is true. plen_float = (self.sf-2)/2 - 2.5*self.has_header + (self.sf-self.ldr*2)/2 * ceil((slen-8)/(self.cr+4)); if no_redundant_bytes plen = ceil( plen_float ); else plen = floor( plen_float ); end end function x_sync = sync(self, x) % sync Packet synchronization % % input: % x: Start index for synchronization % output: % x_sync: Index after up-down alignment % find downchirp found = false; while x < length(self.sig) - self.sample_num up_peak = self.dechirp(x); down_peak = self.dechirp(x, false); if abs(down_peak(1)) > abs(up_peak(1)) % downchirp detected found = true; end x = x + self.sample_num; if found break; end end if ~found return; end % Up-Down Alignment % NOTE preamble_len >= 6 % NOTE there are two NETID chirps between preamble and SFD % NOTE `detect` has already shifted the up peak to position 0 pkd = self.dechirp(x, false); if pkd(2) > self.bin_num / 2 to = round((pkd(2)-1-self.bin_num)/self.zero_padding_ratio); else to = round((pkd(2)-1)/self.zero_padding_ratio); end x = x + to; % set preamble reference bin for CFO elimination pku = self.dechirp(x - 4*self.sample_num); self.preamble_bin = pku(2); if self.preamble_bin > self.bin_num / 2 self.cfo = (self.preamble_bin-self.bin_num-1)*self.bw/self.bin_num; else self.cfo = (self.preamble_bin-1)*self.bw/self.bin_num; end % set x to the start of data symbols pku = self.dechirp(x-self.sample_num); pkd = self.dechirp(x-self.sample_num, false); if abs(pku(1)) > abs(pkd(1)) % current symbol is the first downchirp x_sync = x + round(2.25*self.sample_num); else % current symbol is the second downchirp x_sync = x + round(1.25*self.sample_num); end end function [data_m, checksum_m] = decode(self, symbols_m) % decode Decode data from symbols % % input: % symbols_m: A matrix of symbols to be decoded. Each column % vector represents demodulated symbols from a % LoRa packet. % output: % data_m: A matrix of bytes representing the decoding % result of `symbols_m`. The last two bytes are the % decoded CRC checksum if CRC is enabled. % checksum_m: A vector of checksum based on the decoded % payload.Checksum is empty if CRC is disabled. data_m = []; checksum_m = []; for pkt_num = 1:size(symbols_m, 2) % gray coding symbols_g = self.gray_coding(symbols_m(:, pkt_num)); % deinterleave codewords = self.diag_deinterleave(symbols_g(1:8), self.sf-2); if ~self.has_header nibbles = self.hamming_decode(codewords, 8); else % parse header nibbles = self.hamming_decode(codewords, 8); self.payload_len = double(nibbles(1)*16 + nibbles(2)); self.crc = double(bitand(nibbles(3), 1)); self.cr = double(bitshift(nibbles(3), -1)); % we only calculate header checksum on the first three nibbles % the valid header checksum is considered to be 5 bits % other 3 bits require further reverse engineering header_checksum = [bitand(nibbles(4), 1); de2bi(nibbles(5), 4, 'left-msb')']; header_checksum_calc = self.header_checksum_matrix * gf(reshape(de2bi(nibbles(1:3), 4, 'left-msb')', [], 1)); if any(header_checksum ~= header_checksum_calc) error('Invalid header checksum!'); end nibbles = nibbles(6:end); end rdd = self.cr + 4; for ii = 9:rdd:length(symbols_g)-rdd+1 codewords = self.diag_deinterleave(symbols_g(ii:ii+rdd-1), self.sf-2*self.ldr); % hamming decode nibbles = [nibbles; self.hamming_decode(codewords, rdd)]; end % combine nibbles to bytes bytes = uint8(zeros(min(255, floor(length(nibbles)/2)), 1)); for ii = 1:length(bytes) bytes(ii) = bitor(uint8(nibbles(2*ii-1)), 16*uint8(nibbles(2*ii))); end % dewhitening len = self.payload_len; if self.crc % The last 2 bytes are CRC16 checkcum data = [self.dewhiten(bytes(1:len)); bytes(len+1:len+2)]; % Calculate CRC checksum checksum = self.calc_crc(data(1:len)); else data = self.dewhiten(bytes(1:len)); checksum = []; end data_m = [data_m data]; checksum_m = [checksum_m checksum]; end end function symbols = dynamic_compensation(self, data) % dynamic_compensation Compensate bin drift % % input: % data: Symbols with bin drift % output: % symbols: Symbols after bin calibration % compensate the bin drift caused by Sampling Frequency Offset (SFO) sfo_drift = (1 + (1:length(data))') * 2^self.sf * self.cfo / self.rf_freq; symbols = mod(data - sfo_drift, 2^self.sf); if self.ldr bin_offset = 0; v_last = 1; for i = 1:length(symbols) v = symbols(i); bin_delta = mod(v-v_last, 4); if bin_delta < 2 bin_offset = bin_offset - bin_delta; else bin_offset = bin_offset - bin_delta + 4; end v_last = v; symbols(i) = mod(v+bin_offset, 2^self.sf); end end end function symbols = gray_coding(self, din) % gray_coding Gray coding % `gray_coding` is used in the DECODING process % % input: % data: Symbols with bin drift % output: % symbols: Symbols after bin calibration din(1:8) = floor(din(1:8)/4); if self.ldr din(9:end) = floor(din(9:end)/4); else din(9:end) = mod(din(9:end)-1, 2^self.sf); end s = uint16(din); symbols = bitxor(s, bitshift(s, -1)); self.print_bin("Gray Coding", symbols, self.sf); end function codewords = diag_deinterleave(self, symbols, ppm) % diag_deinterleave Diagonal deinterleaving % % input: % symbols: Symbols after gray coding % ppm: Size with parity bits % output: % codewords: Codewords after deinterleaving b = de2bi(symbols, double(ppm), 'left-msb'); codewords = flipud(bi2de(cell2mat(arrayfun(@(x) ... circshift(b(x,:), [1 1-x]), (1:length(symbols))', 'un', 0))')); self.print_bin("Deinterleave", codewords); end function bytes_w = dewhiten(self, bytes) % dewhiten Data Dewhitening % % input: % bytes: Bytes after deinterleaving % output: % bytes_w: Bytes after dewhitening len = length(bytes); bytes_w = bitxor(uint8(bytes(1:len)), self.whitening_seq(1:len)); self.print_bin("Dewhiten", bytes_w); end function nibbles = hamming_decode(self, codewords, rdd) % hamming_decode Hamming Decoding % % input: % codewords: Codewords after deinterleaving % rdd: Bits with redundancy % output: % nibbles: Nibbles after hamming decoding p1 = LoRaPHY.bit_reduce(@bitxor, codewords, [8 4 3 1]); p2 = LoRaPHY.bit_reduce(@bitxor, codewords, [7 4 2 1]); p3 = LoRaPHY.bit_reduce(@bitxor, codewords, [5 3 2 1]); p4 = LoRaPHY.bit_reduce(@bitxor, codewords, [5 4 3 2 1]); p5 = LoRaPHY.bit_reduce(@bitxor, codewords, [6 4 3 2]); function pf = parity_fix(p) switch p case 3 % 011 wrong b3 pf = 4; case 5 % 101 wrong b4 pf = 8; case 6 % 110 wrong b1 pf = 1; case 7 % 111 wrong b2 pf = 2; otherwise pf = 0; end end if self.hamming_decoding_en switch rdd % TODO report parity error case {5, 6} nibbles = mod(codewords, 16); case {7, 8} parity = p2*4+p3*2+p5; pf = arrayfun(@parity_fix, parity); codewords = bitxor(codewords, uint16(pf)); nibbles = mod(codewords, 16); otherwise % THIS CASE SHOULD NOT HAPPEN error('Invalid Code Rate!'); end else nibbles = mod(codewords, 16); end self.print_bin("Hamming Decode", codewords); end function bytes = symbols_to_bytes(self, symbols) % symbols_to_bytes Derive the payload for a given symbol series % % input: % symbols: Symbols to appear in the demodulation level % output: % bytes: Payload bytes symbols = reshape(symbols, [length(symbols), 1]); self.init(); self.hamming_decoding_en = false; payload_len_ = self.payload_len; if length(symbols) <= 4 slen_tmp = 8 + self.has_header*(self.cr+4); else slen_tmp = 8 + ceil((length(symbols)-4*(1-self.has_header))/4) * (self.cr+4); end self.payload_len = self.calc_payload_len(slen_tmp, true); symbols_ = zeros(self.calc_sym_num(self.payload_len), 1); if self.has_header jj = 9; else jj = 1; end for ii = 1:4:length(symbols) if ii+3 <= length(symbols) symbols_(jj:jj+3) = symbols(ii:ii+3); else symbols_(jj:jj+3) = [symbols(ii:end); zeros(ii-length(symbols)+3, 1)]; end if jj == 1 jj = 9; else jj = jj + self.cr + 4; end end if self.has_header % construct header symbols_tmp = self.encode(zeros(self.payload_len, 1)); symbols_(1:8) = symbols_tmp(1:8); end [bytes, ~] = self.decode(symbols_); if self.crc bytes = bytes(1:end-2); end self.hamming_decoding_en = true; self.payload_len = payload_len_; end function time_ms = time_on_air(self, plen) % time_on_air Calculate the flying time of a LoRa packet % % input: % plen: Payload length % output: % time_ms: Flying time (in milliseconds) sym_num = self.calc_sym_num(plen); % milliseconds time_ms = (sym_num + 4.25 + self.preamble_len) * (2^self.sf/self.bw) * 1000; end function print_bin(self, flag, vec, size) if self.is_debug if nargin == 3 size = 8; end len = length(vec); fprintf("%s:\n", flag); for i = 1:len fprintf("%s\n", dec2bin(round(vec(i)), size)); end fprintf("\n"); end end function print_hex(self, flag, vec) if self.is_debug len = length(vec); fprintf("%s: ", flag); for i = 1:len fprintf("%s ", dec2hex(round(vec(i)))); end fprintf("\n"); end end function log(self, flag, data) if self.is_debug fprintf("%s: ", flag); len = length(data); for i = 1:len fprintf("%d ", data(i)); end fprintf("\n"); end end function plot_peak(self, x) figure; c = [self.downchirp self.upchirp]; for jj = 1:9 for ii = 1:2 ft = fft(self.sig(x:x+self.sample_num-1).*c(:,ii), self.fft_len); ft_ = abs(ft(1:self.bin_num)) + abs(ft(self.fft_len-self.bin_num+1:self.fft_len)); subplot(2, 9, (ii-1)*9+jj); plot(ft_); end x = x + self.sample_num; end end end methods(Static) function b = bit_reduce(fn, w, pos) b = bitget(w, pos(1)); for i = 2:length(pos) b = fn(b, bitget(w, pos(i))); end end function w = word_reduce(fn, ws) w = ws(1); for i = 2:length(ws) w = fn(w, ws(i)); end end function y = topn(pks, n, padding, th) [y, p] = sort(abs(pks(:,1)), 'descend'); if nargin == 1 return; end nn = min(n, size(pks, 1)); if nargin >= 3 && padding y = [pks(p(1:nn), :); zeros(n-nn, size(pks, 2))]; else y = pks(p(1:nn), :); end if nargin == 4 ii = 1; while ii <= size(y,1) if abs(y(ii,1)) < th break; end ii = ii + 1; end y = y(1:ii-1, :); end end function y = chirp(is_up, sf, bw, fs, h, cfo, tdelta, tscale) % chirp Generate a LoRa chirp symbol % % input: % is_up: `true` if constructing an up-chirp % `false` if constructing a down-chirp % sf: Spreading Factor % bw: Bandwidth % fs: Sampling Frequency % h: Start frequency offset (0 to 2^SF-1) % cfo: Carrier Frequency Offset % tdelta: Time offset (0 to 1/fs) % tscale: Scaling the sampling frequency % output: % y: Generated LoRa symbol if nargin < 8 tscale = 1; end if nargin < 7 tdelta = 0; end if nargin < 6 cfo = 0; end N = 2^sf; T = N/bw; samp_per_sym = round(fs/bw*N); h_orig = h; h = round(h); cfo = cfo + (h_orig - h) / N * bw; if is_up k = bw/T; f0 = -bw/2+cfo; else k = -bw/T; f0 = bw/2+cfo; end % retain last element to calculate phase t = (0:samp_per_sym*(N-h)/N)/fs*tscale + tdelta; snum = length(t); c1 = exp(1j*2*pi*(t.*(f0+k*T*h/N+0.5*k*t))); if snum == 0 phi = 0; else phi = angle(c1(snum)); end t = (0:samp_per_sym*h/N-1)/fs + tdelta; c2 = exp(1j*(phi + 2*pi*(t.*(f0+0.5*k*t)))); y = cat(2, c1(1:snum-1), c2).'; end function spec(sig, fs, bw, sf) x = 0:1/sf:2^sf/bw; y = -bw/2:bw/2; p = fs/bw; win_length = 2^(sf-2); N = p*2^sf; s = spectrogram(sig,win_length,round(win_length*0.8),N); valid_data_len = round(2^sf/2*1.5); b = abs(s(valid_data_len:-1:1,:)); c = abs(s(end:-1:end-valid_data_len+1,:)); d = [b;c]; figure; imagesc(x,y,d); colormap summer; title('Spectrogram'); xlabel('Time'); ylabel('Frequency'); set(gcf,'unit','normalized','position',[0.05,0.2,0.9,0.1]); end % Method `read` and `write` are copied from % https://github.com/gnuradio/gnuradio/blob/master/gr-utils/octave/read_complex_binary.m % https://github.com/gnuradio/gnuradio/blob/master/gr-utils/octave/write_complex_binary.m function v = read(filename, count) % usage: read(filename, [count]) % % open filename and return the contents as a column vector, % treating them as 32 bit complex numbers % m = nargchk (1,2,nargin); if (m) usage (m); end if (nargin < 2) count = Inf; end f = fopen (filename, 'rb'); if (f < 0) v = 0; else t = fread (f, [2, count], 'float'); fclose (f); v = t(1,:) + t(2,:)*1i; [r, c] = size (v); v = reshape (v, c, r); end end function v = write(data, filename) % usage: write(data, filename) % % open filename and write data to it % Format is interleaved float IQ e.g. each % I,Q 32-bit float IQIQIQ.... % This is compatible with read_complex_binary() % m = nargchk (2,2,nargin); if (m) usage (m); end f = fopen (filename, 'wb'); if (f < 0) v = 0; else re = real(data); im = imag(data); re = re(:)'; im = im(:)'; y = [re;im]; y = y(:); v = fwrite (f, y, 'float'); fclose (f); end end end end