USRP变速率采样原理及实现

公司网址：www.luowave.com   国内首发   USRP LW-X310

USRP是数款流行的SDR硬件中功能和应用都相对成熟的一款产品，从WIFI协议、ZigBee协议、RFID协议、GSM通信系统、LTE 4G通信系统到飞机通信、卫星通信USRP都能很好的进行支持。软件开发工程师可以用它开发应用，安全工程师则用它来测试、研究相关的无线通信协议

1. 变速率采样

5M带宽的LTE系统中，采样率为512＊15k＝7.68MHz，而USRP2在不改变FPGA代码的情况下仅支持100M的整数倍分频的采样率，最接近7.68M的采样率为6.25M。所以在空口获得的数据采样率为6.25M，而LTE系统内部采样率为7.68M，需要进行插值和抽取实现速率转换。

（补充说明一点，由于LTE过采样的特点，采样率6.25M仍然可以完整保留信号特性，实际上，只要大于5M的采样率都是可以的，7.68M应该是为了作FFT的方便设定的）

以接收为例，如果进行6.25M到7.68M的采样率转换，需要先进行768倍内插，再进行625倍抽取，即6.25M * 768 / 625 = 7.68M，这样可以防止丢失采样点信息（先抽取后内插会丢失）。通用的变速率流程是：

插767个0（x00...0 x00..0 x...）——LPF低通滤波——抽取(每625个点抽取一个)

（插0后过LPF实际上完成了插值的功能，而在抽取前过LPF则是为了防止抽取带来的镜像干扰，两个LPF可以合并为一个，截止频率用低的即可）

2. 低通滤波器设计

低通滤波的作用是消除内插和抽取带来的镜像干扰，低通截止频率为min(7.68,6.25)/2=3.125MHz，然而对于5M带宽的信号来讲，基带信号带宽为2.5MHz，实际上在2.5~3.125之间并没有信号分量，所以在设计滤波器的时候可以放宽过渡带带宽的要求，通带范围是0～2.5Hz，截止频率即阻带起始频率为3.125+3.125-2.5=3.75MHz。阻带衰减根据需要而定，默认80db。注意滤波器是在插值之后，故滤波器的Fs＝6.25M*768=4800MHz，这样设计的滤波器阶数一般很高，程序中采用10000阶，据说超过10000阶fdatool会跑死，有兴趣可以试验一下。

 

利用matlab fdatool可以进行滤波器设计，浮点系数到定点系数的转换。在matlab命令行中输入”fdatool”即可打开设计面板，参数介绍参考文档《MATLAB FDATOOL的使用.doc》。

（补充一点，对于定点和浮点的转换主要考虑到计算效率的问题，设计定点滤波器的时候需要考虑溢出和精度的问题）

 

分析fir滤波器幅频特性和相频特性

b=[fir系数向量]

a=1

fvtool(b,a)

或者freqz(b,a,256,fs)

 

上述操作均在windows下matlab完成，octave没试过。

3. 变速率采样的实现——多相滤波器

上面给出的变速率采样通用流程中，插值，滤波，抽取，这个过程实现复杂度很高，以接收为例，一个subframe的空口数据有6250个sample，插0后变成6250＊768＝4800000个sample，滤波器抽头数目为10000，滤波的计算量非常大，实际中采用多相滤波器来简化运算，对于多相分解的原理不用深究，其思想主要有两点：

(1)由于抽取中需要舍弃多个点，例如每625个点抽取一个，对于没抽到的点，不需要进行卷积运算

(2)滤波器的输入向量中含有大量的0，这些0不会影响输出结果，其实只有插值前的数据会与某几个抽头进行乘法累加的运算，卷积操作可大大简化。

根据上面分析，多相滤波器的关键在于，对于每一个输出采样点（抽取后的），找到其对应的输入相位，

即输入向量的偏移位置，再对其中不为0的样点进行乘法累加运算。

 

实现程序如下，P为升采样倍数768，Q为降采样倍数625（发送相反）：

 

第一步，计算出输出采样点的数目，初始化输入向量和输出向量以及滤波器向量的起止位置。

所有向量都是按照IQIQIQ...的形式，I和Q都是short型。

    out_sample_count = (int)(in_sample_count * (float)P / (float)Q);

    short *output_ptr = out_data;

    unsigned int output_ix = 0;

    short const *input_end = in_data + in_sample_count * 2;

    short const *output_end = out_data + out_sample_count * 2;

    short const *filter_end = filter_coeff + filter_len;

 

第二步，对于每一个输出采样点，进行乘法累加运算（卷积实现滤波）。其中关键步骤用abcd标出。

    while (output_ptr < output_end) {

 

2－a. 对于特定输出采样点，确定输入向量偏移位置以及对应的滤波器抽头下标。为了省去这部分计算量，采用查表的形式，table在初始化阶段建立。

        //branch_table[i] = (i * Q) % P;

        //offset_table[i] = (i * Q - branch_table[i]) / P;

 

        int output_branch = branch_table[output_ix];

        int input_offset = offset_table[output_ix];

       

        short *input_ptr = in_data + 2 * input_offset;

        short *filter_ptr = filter_coeff + output_branch;

       

        while (input_ptr >= input_end) {

            input_ptr -= 2;

            filter_ptr += P;

        }

 

2－b. 卷积运算关键，乘法累加，input_ptr -= 2是因为输入向量在卷积时要先做反转再做乘法累加，而filter_ptr += P则是因为其中忽略了插值0对应的抽头。

       

        int sum_real = 0, sum_imag = 0;

        while ((input_ptr >= in_data) && (filter_ptr < filter_end)) {

            sum_real += (int)(*input_ptr) * (int)(*filter_ptr);

            sum_imag += (int)(*(input_ptr + 1)) * (int)(*filter_ptr);

            input_ptr -= 2;

            filter_ptr += P;

        }

 

2－c. 为防止溢出，sum的IQ都是int型的，在输出向量为short的时候需要进行缩放，即移位操作

        *output_ptr = (short)(sum_real >> fix_shift);

        *(output_ptr + 1) = (short)(sum_imag >> fix_shift);

 

        output_ptr += 2;

        output_ix++;

    }

举个例子，假如我们计算输出的第5个采样点的数据，那么实际上抽取前的输出采样点下标应该是625＊4＝2500（output_idx=4，假设向量下标都是从0开始，抽取的第一个采样点是output[0]），这个点对应的计算情况

      h0 …h196，....，h964，....，h1732,....，h2500

0..0，0..0，x3，767个0，x2，767个0，x1，767个0，x0

 572  195   1   ｜          768＊3＝2304个sample            ｜

      ＊ 2304＋196＝2500，＊即为第2500个点对应的输入向量位置

 

如果进行乘法累加，实际上只需要计算四个点即可，首先找到最近的非零输入点为x3（［4＊625－196］／768），滤波器抽头h196（4＊625％768），而后只需要计算非零输入点对应的乘法累加即可。在代码里体现为：

            sum_real += (int)(*input_ptr) * (int)(*filter_ptr);

            sum_imag += (int)(*(input_ptr + 1)) * (int)(*filter_ptr);

            input_ptr -= 2;

            filter_ptr += P;

 

这样既不需要插0，也不需要对舍弃点进行滤波运算。

 

4. 无限长序列滤波——重叠保留法实现的分段卷积

因为输入序列是从LTE协议栈产生的（TX）或者从空口接收到的（RX），所以序列都是无限长的，所以在滤波时实质上是用的分段卷积。保留的历史数据长度只要满足L*P>滤波器长度，滤波后再舍弃这部分历史数据产生的结果就可以啦。

以发送为例

4-a. 历史数据

memcpy(inData, inHistory, IN_HISTORY_LEN*4);

4-b. 追加上下一段序列

memcpy(inData+IN_HISTORY_LEN*2, inSignal, IN_SIGNAL_LEN*4);

4-c. 滤波

polyphaseResample(inData, IN_DATA_LEN, outData, OUT_DATA_LEN, coeff_3072_to_25, sizeof(coeff_3072_to_25) / sizeof(short), OUTRATE, INRATE, 11, g_branch_table, g_offset_table);

4-d. 舍弃历史数据产生的输出，返回滤波结果

memcpy(gTxResampledDataVector+iSubFrameIdx*OUT_SIGNAL_LEN*2, outData+OUT_HISTORY_LEN*2, OUT_SIGNAL_LEN*4);

4-e. 保留部分输入序列作为历史数据

memcpy(inHistory, inSignal+IN_SIGNAL_LEN*2-IN_HISTORY_LEN*2, IN_HISTORY_LEN*4);

4 仿真结果测试

原始信号波形

变采样信号（7.68M->6,25M）

会发现过滤波器有一定时延，不过这个时延是针对发送信号做得整体的偏移，不会影响后面的发送数据。

对于原始信号，先进行7.68->6.25M变采样，然后进行6.25->7.68M变采样，利用同步模块解调，证明变速率的过程是可逆的，不会丢失信号信息。