在上一篇文章中,已经完成了BPSK波形的发射。
相对于BPSK波形的生成总共就4行代码,接收要略微复杂一些,算上各种同步、锁相环,约80行。完整版参考Git仓库。
设备连接:
USRP B210 基于 AD9361芯片,把放大、变频、滤波、AD、数字变频都做完了。输入的是射频信号,输出的就是基带(0中频)复平面上的X、Y坐标。因此,解调代码输入的就是XY坐标(IQ路信号)。这里给出其C语言核心代码如下:
vector<char> demod_psk(const short (*iqdata)[2], const int splen)
{static const double fir_rcos_q25[25] = {-0.018773,0.0030136,0.032677,0.047094,0.02655,-0.027522,-0.085225,-0.099447,-0.032147,0.11904,0.31118,0.472,0.53463,0.472,0.31118,0.11904,-0.032147,-0.099447,-0.085225,-0.027522,0.02655,0.047094,0.032677,0.0030136,-0.018773};static double fir_cache_i[25] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};static double fir_cache_q[25] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};//fir filterstatic unsigned long long clk = 0, bad_angle = 0, last_synclk = 0;static const int WINSZ = 4000;static double stati_win[WINSZ][2], sum_win[2] = {0,0};static double cache_win[WINSZ][2];//smp best posstatic double clk_offset[4] = {0,0,0,0};static bool bInited = false;static double curr_dp = 0;static int bestOff = 0;//stati_win,cache_win 归零初始化略去//...//bpsk bit outvector<char> dp_bitbuf;//逐个点步进for (int i=0;i<splen;++i){//1. 滤波去带外噪声fir_cache_i[clk % 25] = iqdata[i][0];fir_cache_q[clk % 25] = iqdata[i][1];double fval_i = 0,fval_q = 0;for (int f=0;f<25;++f){fval_i += fir_cache_i[(clk+f+1) % 25] * fir_rcos_q25[f];fval_q += fir_cache_q[(clk+f+1) % 25] * fir_rcos_q25[f];}//2. 窗口统计去直流分量sum_win[0] -= stati_win[clk%WINSZ][0];sum_win[1] -= stati_win[clk%WINSZ][1];sum_win[0] += fval_i;sum_win[1] += fval_q;stati_win[clk%WINSZ][0] = fval_i;stati_win[clk%WINSZ][1] = fval_q;const double dcremoved[2] {fval_i - sum_win[0]/WINSZ, fval_q - sum_win[1]/WINSZ };cache_win[clk%WINSZ][0] = dcremoved[0];cache_win[clk%WINSZ][1] = dcremoved[1];//3. 最佳采样位置跟踪const int curr_offv = clk % 4;clk_offset[curr_offv] += (dcremoved[0]*dcremoved[0] + dcremoved[1] * dcremoved[1]); //符号时钟位置同步double phy = 0;//锁相环误差if (curr_offv == bestOff){//锁相旋转double rt[2]{cos(-curr_dp),sin(-curr_dp)};const double out[2]{cache_win[(clk+4)%WINSZ][0],cache_win[(clk+4)%WINSZ][1]};double rp[2]{out[0] * rt[0] - out[1] * rt[1],out[0] * rt[1] + out[1] * rt[0]};//失锁重置,即重新测量相位误差phy = atan_r(rp[1],rp[0]);if (phy >= c_pi/4 || phy <=-c_pi/4){curr_dp = atan_r(out[1],out[0]);rt[0] = cos(-curr_dp);rt[1] = -sin(curr_dp);rp[0] = out[0] * rt[0] - out[1] * rt[1];rp[1] = out[0] * rt[1] + out[1] * rt[0];phy = atan_r(rp[1],rp[0]);++bad_angle;}//更新环路curr_dp += (phy/2);curr_dp = wrapto2pi(curr_dp );//解调后的01比特dp_bitbuf.push_back( rp[0] < 0? 0: 1 );}//错开最佳点2符号进行时钟误差补偿if (last_synclk + WINSZ < clk && ((clk+2) % 4)==bestOff){last_synclk = clk;int maxo = 0;double maxv = clk_offset[0];if (clk_offset[1] > maxv) { maxo = 1; maxv = clk_offset[1];}if (clk_offset[2] > maxv) { maxo = 2; maxv = clk_offset[2];}if (clk_offset[3] > maxv) { maxo = 3; maxv = clk_offset[3];}bestOff = maxo;clk_offset[0] = 0; clk_offset[1] = 0;clk_offset[2] = 0; clk_offset[3] = 0;}++clk;}return dp_bitbuf;
}
不停地调用上述代码,输入一段一段接收到的X\Y坐标(IQ路信号),输出一串串01比特。上述代码采用全静态内存、无任何的搬移、拷贝,基本满足一般低速信号的处理要求。
下面我们逐一介绍每一步干了什么。在此之前,先讨论一个题外话。
0. 题外话:软硬件的思维大不相同
先放一张图:
有多少人在看到这个BPSK 科斯塔斯环时,会思考这个东西的C语言如何实现VCO、LF、LPF,这么多滤波如何提速等等。其实我也困扰了很久.,其实若从复平面、X\Y坐标去理解,会发现用PC计算机实现BPSK解调,完全不必有VCO以及这些设施。
上面的示意图,是在上个世纪使用模拟器件处理解调工作的最佳方案。硬件记不住太多的历史数据,因此无法统计一些信息;硬件无法方便的计算各种角度,比如atan。所以它用一个压控震荡源产生一个对消载波,去环路匹配残差频率。在C语言里,我们可以直接测量角度,从而跟踪上向量的旋转,对消载波是不必要的。PC的思路是:
符号时刻t0,观察到向量(x,y)在131度,对BPSK,则下一个符号时刻 t1 要么是 131+180,要么还是131度。因此,直接测量 atan(y,x),而后稍微和预计位置比较一下,就可以得到残差相位了,太直接了。整个环路就等效为一个角度变量的跟踪。
旋转打地鼠
如果把接收BPSK的两个可能的符号位置想象成一个黑色幕布下的两组地鼠洞,地鼠会随时从任意洞口钻出来。更要命的是,整个地鼠机放在一个转台上,转台的转速就是各级模拟/数字下变频(ADC、DDC)到0频(基带,或者复平面xy坐标)时保留的频率误差。现在,要准确的打地鼠,就是和解调非常类似的活动。
按照4倍速率XY坐标运动的地鼠 | 角度跟踪最佳状态的地鼠 |
---|---|
传统教材的科斯塔斯环,等于是要产生一个和残差相位相反的旋转矢量,叠加到输入载波,让整体的地鼠盘子不动。而我们的C语言解调,只要测量当前的残差相位,并不断跟踪即可。盘子继续转动,多了一个测量杆子(curr_dp),这个杆子和地鼠一起转动,每次根据地鼠的位置,跟踪调整杆子的角度即可。
1. 滤波去带外噪声
由于要进行时钟同步,以及奈奎斯特限制,接收时往往采用比发射速率高几倍的速率进行接收。此时,信号带宽外的噪声也进来了。去除带外噪声,可以一定程度良化波形。
上文的代码:
//1. 滤波去带外噪声fir_cache_i[clk % 25] = iqdata[i][0];fir_cache_q[clk % 25] = iqdata[i][1];double fval_i = 0,fval_q = 0;for (int f=0;f<25;++f){fval_i += fir_cache_i[(clk+f+1) % 25] * fir_rcos_q25[f];fval_q += fir_cache_q[(clk+f+1) % 25] * fir_rcos_q25[f];}
构造了一个环形结构,实现时域滤波的卷积计算。由于FIR是对称的,导致卷积反转不再需要。只要用时钟咬紧当前的位置即可,搞一个25长度的环状缓存,只记忆滤波需要的部分。同时,时钟推进后,用取余替代memcpy,不进行内存搬移,几乎没有开销。
这部分的效果,可以通过等效的Octave代码处理真实数据来展示:
左图是滤波前的向量位置打点图,右图是滤波后的。显然滤波后,0、1两坨坐标分的更开了。下面是仿真代码:
pkg load communications;
clear;
clc;
mulrate = 4;
symlen = 1000;
dta1_bits = randint (1, symlen, 2);
dta2_symbols = pskmod (dta1_bits, 2, 0, "gray");
#dta2_symbols = qammod (dta1_bits, 4);
#figure(); plot (dta2_symbols, "bx");dta3_sig = zeros(1,symlen*mulrate);
dta3_sig(1:mulrate:end) = dta2_symbols;#figure(); plot (dta3_sig, "r+");##fir send
[fir_rcos]=rcosfir(0.25,[-3,3],mulrate,1,'sqrt');dta4_baseband = conv(dta3_sig,fir_rcos,'same');#plot (dta4_baseband(1:mulrate:end), "g*");
#figure; plot(fir_rcos);
##Trans# to rx
rf_sym_rate = 1000;
rf_sample_rate = rf_sym_rate * mulrate;dta5_recv = awgn (dta4_baseband, 12);figure();plot (dta5_recv(1:mulrate:end), "r.");##fir recv
dta6_filted = dta5_recv;
dta6_filted = conv(dta5_recv,fir_rcos,'same');figure();plot (dta6_filted(1:mulrate:end), "b.");dta7_sampled = dta6_filted([1:mulrate:end]);
#figure();plot (dta7_sampled, "b.");#dta8_est = qamdemod (dta7_sampled, 4);
dta8_est = pskdemod (dta7_sampled, 2,0,"gray");biterr (dta1_bits, dta8_est)#figure();plot (real(dta4_baseband),'r');
#figure();plot (real(dta6_filted),'b');
2 窗口统计去直流分量
有些时候,收到的向量具有一个整体的直流分量,导致其不再围绕原点分布。表现在复平面上,就是一坨位置跑到别处去。这样导致基于原点进行角度测量的三角函数都没有办法使用了。
解决此问题的方法很简单,就是统计一个窗口里的坐标的均值,并搬移到原点去。上文代码:
//2. 窗口统计去直流分量sum_win[0] -= stati_win[clk%WINSZ][0];sum_win[1] -= stati_win[clk%WINSZ][1];sum_win[0] += fval_i;sum_win[1] += fval_q;stati_win[clk%WINSZ][0] = fval_i;stati_win[clk%WINSZ][1] = fval_q;const double dcremoved[2] {fval_i - sum_win[0]/WINSZ, fval_q - sum_win[1]/WINSZ };
就干了这个事情。等一等,为什么没有for循环?这里面又使用了第二个环状统计缓存器,即均值统计环状缓存。
- 缓存记住最新的WINSZ = 4000个坐标,以及他们的和sum_win。
- 新的数据到来,sum_win 先减去最老的那个时刻的值,并加上新时刻的值。
- 如此循环,每个时钟下,窗口都能跟踪当前的均值了。
3 时钟和频率同步
这是最关键的部分,首先是时钟同步,而后是频率同步。
时钟同步是为了获得地鼠盘,让地鼠稳定的出现在一个看不见的圆圈上,而不是乱窜。
频率同步就是打地鼠,用一个角度跟踪器(软件实现)或者锁相环(硬件)跟踪残差频率造成的相位移动。
3.1 在恰当的时刻获取符号
在4倍采样率下,由于不知道下图红色的正确位置在哪里,所以随便抽取一路作为x,y的坐标值是不行的。
如何获取正确的位置呢?我们通过观察,这些正确位置上,矢量均远离坐标原点;在过渡位置上,一半以上的点贴近坐标原点。因此,只要统计当前四个偏移位置(01,2,3)上,4000个样点的矢量长度和,和最大的哪个就是最佳位置。
//3. 最佳采样位置跟踪const int curr_offv = clk % 4;clk_offset[curr_offv] += (dcremoved[0]*dcremoved[0] + dcremoved[1] * dcremoved[1]); //符号时钟位置同步
这里求矢量长度没有取SQRT,为了省计算。其实是一样的,只是比大小,sqrt不改变其单调性。
由于收发双方的时钟存在误差,这个最佳位置是会不断移动的。一小会是2,后面可能就成了3,0,1, 或者反方向编程了1,0,3,2的重复。所以,每隔一段时间,我们要重新评估这个最佳位置,并只在最佳位置上抽取。其逻辑片段:
static unsigned long long clk = 0, last_synclk = 0;static int bestOff = 0;//逐个点步进for (int i=0;i<splen;++i){//...//3. 最佳采样位置跟踪const int curr_offv = clk % 4;clk_offset[curr_offv] += (dcremoved[0]*dcremoved[0] + dcremoved[1] * dcremoved[1]); //符号时钟位置同步//在最佳位置输出if (curr_offv == bestOff){//...//解调后的01比特dp_bitbuf.push_back( rp[0] < 0? 0: 1 );}//错开最佳点2符号进行重新评估if (last_synclk + WINSZ < clk && ((clk+2) % 4)==bestOff){last_synclk = clk;int maxo = 0;double maxv = clk_offset[0];if (clk_offset[1] > maxv) { maxo = 1; maxv = clk_offset[1];}if (clk_offset[2] > maxv) { maxo = 2; maxv = clk_offset[2];}if (clk_offset[3] > maxv) { maxo = 3; maxv = clk_offset[3];}bestOff = maxo;clk_offset[0] = 0; clk_offset[1] = 0;clk_offset[2] = 0; clk_offset[3] = 0;}++clk;}
错开最佳点2符号进行重新评估,是为了稳定、隐含的解决双方钟差不同步需要进行符号填补、扣除的操作。比如收方每收10000个符号,发方发了10017或者9983个符号,那就要悄默默的追赶/剔除17个,否则就乱套了。这种追赶/剔除行为,就是在 bestOff调整时顺带完成的。bestOff变了,导致下一次输出 if (curr_offv == bestOff) 判断会提前、退后一个点。为了这个行为的稳定,对bestOff的调整要避开当前最佳点,越远越好,于是时机选择在 ((clk+2) % 4)==bestOff。
3.2 测量并咬住相位
一旦找到了正确的时钟位置,测量相位就是一个 x,y坐标到角度的变化了:
static unsigned long long clk = 0, last_synclk = 0;static int bestOff = 0;static double curr_dp = 0;double phy = 0;//锁相环误差if (curr_offv == bestOff){//锁相旋转double rt[2]{cos(-curr_dp),sin(-curr_dp)};const double out[2]{cache_win[(clk+4)%WINSZ][0],cache_win[(clk+4)%WINSZ][1]};double rp[2]{out[0] * rt[0] - out[1] * rt[1],out[0] * rt[1] + out[1] * rt[0]};//失败重置//...//更新环路curr_dp += (phy/2);curr_dp = wrapto2pi(curr_dp );//解调后的01比特dp_bitbuf.push_back( rp[0] < 0? 0: 1 );}//...++clk;}
当然,如果角度误差巨大,说明跟踪失败了,可以重置:
//失锁重置,即重新测量相位误差phy = atan_r(rp[1],rp[0]);if (phy >= c_pi/4 || phy <=-c_pi/4){curr_dp = atan_r(out[1],out[0]);rt[0] = cos(-curr_dp);rt[1] = -sin(curr_dp);rp[0] = out[0] * rt[0] - out[1] * rt[1];rp[1] = out[0] * rt[1] + out[1] * rt[0];phy = atan_r(rp[1],rp[0]);++bad_angle;}
4. 最终效果
通过上述操作,最终紧紧咬住了BPSK 180度的相位变化,获得了正确的数据流。由于不知道确切的相位是0还是180,可能存在极性颠倒。如何抗颠倒,交给上层协议去做了。
完整的SDR实现了双工的协议栈,用主对方两个频率实现全通的IP网络连接。对于从0开始实现的纯C/C++ SDR大玩具,可以非常细致的体会真正的通信系统设计实现环节的各种坑。目前还没有涉及到精确定时。精确定时在粗大的颗粒度下,SDR是可以来做的。要在20ms以下进行定时、TDMA操作,恐怕需要再踩一些坑。不过能够用有限的知识,控制低成本的实验设备完成全套流程,可玩性已经非常高了。
完整代码参考前文的Git仓库。发行版参考前文连接。
5.思考体会
通过这个实验,我们发现对于基带信号,直接从复平面的特性去思考,要比从受限于模拟器件的角度得到的原理图去思考更直接。使用通用计算机进行操作,或许有更为简便的算法等待开发。应该有教材在介绍传统的解调思路前,把这些思考说清楚,这样学生在学习时,就更为有针对性了。