USRP变速率采样原理及实现

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可以合并为一个，截止频率用低的即可）
低通滤波的作用是消除内插和抽取带来的镜像干扰，低通截止频率为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没试过。
上面给出的变速率采样通用流程中，插值，滤波，抽取，这个过程实现复杂度很高，以接收为例，一个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 * Q) % P;
        //offset_table = (i * Q - branch_table) / 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，也不需要对舍弃点进行滤波运算。

因为输入序列是从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);
原始信号波形

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

会发现过滤波器有一定时延，不过这个时延是针对发送信号做得整体的偏移，不会影响后面的发送数据。
对于原始信号，先进行7.68->6.25M变采样，然后进行6.25->7.68M变采样，利用同步模块解调，证明变速率的过程是可逆的，不会丢失信号信息。