多相抽取滤波原理与定点化实现（Matlab）_matlab 多项滤波抽取

1.多相抽取滤波原理

204B输出四相数据多项滤波抽取过程可表示为：

其中， x(n)与y(n) 分别为输入与输出信号；h(n)为多项滤波器的冲激响应。

输入可以写作

上式表示将输入分别延时1，2…M-1点，并抽取M倍后的信号。

类似地，多相滤波器的z变换可以写作

即：

令每一行的和为：

有

与输入相卷结构如下：

注意到输入没有做抽取，只是在每个支路上依次延时了一个符号。在完成前面的等效滤波部分后，合成一路，再进行抽取。

利用诺贝尔恒等变换进行抽取与延时的交换：

其中

这个结构相当于将输入的1、2、3…M-1整数倍时刻的点依次送入0、1、2…M-1子滤波器后求和。

而输入与输出的支路并不是按顺序一一对应的，这是因为延迟后的抽取会使拿到的数据与我们常识中的顺序不同。

用0来表示输入的延时，可以发现第二相并不是2、5，因为数据延时后抽取到的第一个数在最后一个位置。因此常识中的“第二相“会出现在最后一相。

注意到无论在输入的第几相延时几拍，输出会固定延时一拍；这是因为无论输入数据延迟几拍，在分解为多相数据后相当于每一相都进行了抽取，除了第一相外多相数据都会得到一拍的固定延时。

输入数据不存在延时时，可以选择将输出的除第一相外的数据延时一拍求和得到。其等效性很显然。

关于原理，弗雷德里克·J·哈里斯的《通信系统中的多采样率信号处理》一书有很详细的推导，建议阅读。

2. Matlab验证

        设计目标：512M采样384M中频数据，带宽为100M；128M混频到基带；不过低通，因为下一级需要进行抽取四倍，需要经过一个4相的低通，抽取完成后采样率为128M。

        在matlab进行定点化，将使用/不使用多项结构的抽取滤波结果以及verilog的计算结果对比，要求结果完全一致。

2.1数据产生

首先产生512M采样384M信号，信号为384M左右分别10M、20M的双音信号，采样后定点化为14位有符号数。

这个过程如下，根据带通采样定理可知，384M中心频率的信号被搬到了128M

2.2 正交混频

产生正交的载波将信号频谱平移：

但注意到混频NCO频率为128M，而采样率为512M，因此得到的载波仅有3个取值：0、-1、1。因此这部分计算可以省略，并通过简单的保持、取反、取0得到。这部分并没有精度的变化，因此信号格式依旧为s(14,0)。

2.3多相分解

         将混频数据拆为4相。实际上混频也是多相进行的，本次设计中混频的计算省略了，因此混频前后拆为4相对应的实现结构并没有区别

二三四支路延时一拍，与系数对应关系：1-1；2-4；3-3；4-2。

得到结果与直接滤波抽取完全一致。

3.FPGA实现与验证

FPGA使用以下结构实现，128点系数拆解为4相，每相32点；注意到每个子滤波器的输出先与其它几项求和再进行截断与饱和，这样计算过程才与直接型FIR等效。

结果对比：考虑到滤波器带来的延迟，输出的有效数据在matlab卷积结果的第32位开始保持一致。

进行进一步的对比

得到最大偏差为0，即matlab运算结果与FPGA运算结果保持完全一致。

4.Matlab源码

主文件：

clc;
clear all;
 
%% 参数定义
fc1  = 384e6;      %输入中心频率
fc2  = 128e6;      %混频频率
fs   = 512e6;      %采样率
span = 80e6;      %带宽
t    = 0:512*10-1;      %采样5120点
N    = length(t);
%y_sample = 2^12*((cos(2*pi*(fc1-span/8)/fs.*t)+cos(2*pi*(fc1+span/4)/fs.*t)));
y_sample = 0.5*((cos(2*pi*(fc1-span/8)/fs.*t)+cos(2*pi*(fc1+span/4)/fs.*t)));
%定点化为s[14,0]
qpath  = quantizer('fixed','round','saturate',[14,0]);
fix_y  = quantize(qpath,2^13*y_sample);
err1 = fix_y-2^13*y_sample;
%% 混频
% mix 128MHz sampled by 512MHz 
% nco_sin : 0 1 0 -1
% nco_cos : 1 0 -1 0
%IQ swapped
phase1 = upsample(fix_y(1:4:end),4);
phase2 = upsample(fix_y(2:4:end),4);
phase3 = upsample(fix_y(3:4:end),4);
phase4 = upsample(fix_y(4:4:end),4);
mix_q = phase1-[0 0 phase3(1:end-2)];
mix_i = -[0 phase2(1:end-1)]+[0 0 0 phase4(1:end-3)];
mix_out = mix_i + 1j* mix_q;
 
nco_cos = 2^15*cos(2*pi*(fc2)/fs.*t);
nco_sin = 2^15*cos(2*pi*(fc2)/fs.*t-pi/2);
qpath2  = quantizer('fixed','round','saturate',[16,0]);
fix_nco_cos  = quantize(qpath2,nco_cos);%fix_y s[14,0]
fix_nco_sin  = quantize(qpath2,nco_sin);%fix_y s[14,0]
%混频
fix_y_q = fix_nco_cos.*fix_y;%s[14,0] * s[16,0] s[29,0]
fix_y_i = -fix_nco_sin.*fix_y;
fix_y_nco_out = fix_y_i+1j*fix_y_q;
fix_y_nco_out = round(fix_y_nco_out./2^15);%得到s[14,0]
%拆解为4相数据
              mix_out_phase1=mix_out(1:4:end);
              mix_out_phase2=mix_out(2:4:end);
              mix_out_phase3=mix_out(3:4:end);
              mix_out_phase4=mix_out(4:4:end);
 
path = 'D:\proj\matlab_proj\DDC\DDC\IQ_mixer+Poly+CIC+HB+FIR Span=100MHz\PolyDecim0_Coef_Fp0.19531.coe';
Poly0 = coeread(path);
Poly0_coe = Poly0.Numerator;
Poly0_coe_p1 = Poly0_coe(1:4:end);    
Poly0_coe_p2 = Poly0_coe(2:4:end);    
Poly0_coe_p3 = Poly0_coe(3:4:end);    
Poly0_coe_p4 = Poly0_coe(4:4:end);    
%% 四路 分别滤波并抽取
y_Poly1=[conv(mix_out_phase1,Poly0_coe_p1) 0];%补个0好加起来
y_Poly2=[0 conv(mix_out_phase2,Poly0_coe_p4)];
y_Poly3=[0 conv(mix_out_phase3,Poly0_coe_p3)];
y_Poly4=[0 conv(mix_out_phase4,Poly0_coe_p2)];
%将四路加起来
y_Poly5=y_Poly1+y_Poly2+y_Poly3+y_Poly4;%4*s[35,15]=s[37,15]
y_Poly5_rnd = round(y_Poly5./2^15); %去除小数 s[22,0]
y_Poly5_sat = sat(y_Poly5_rnd,16,1);  %去除6位整数 s[16,0]
%直接滤波抽取，不做多相分解
%% 直接滤波抽取，不做多相分解
y_Poly0=FIR_deci(mix_out,Poly0_coe,4,15);
max(abs(y_Poly5_sat-y_Poly0))%说明两种计算等价
 
%% FPGA 结果对比
y_Poly0_i = real(y_Poly0);
y_Poly0_q = imag(y_Poly0);
 
y_Poly0_i = y_Poly0_i(32:end);
y_Poly0_q = y_Poly0_q(32:end);
poly_i = load("2048_i.txt")';
poly_q = load("2048_q.txt")';
 
err_i = poly_i - y_Poly0_i(1:length(poly_i));
err_q = poly_q - y_Poly0_q(1:length(poly_q));
max(abs(err_i))
max(abs(err_q))
 
%% 滤波结果
figure(1)
subplot(5,1,1)
plot_fft(Poly0_coe,512,1,1); title('滤波器响应')
subplot(5,1,2)
plot_fft(y_Poly5_sat,128,1,1);title('多相滤波后频谱')
subplot(5,1,3)
plot(real(y_Poly5_sat));title('多相滤波后数据')
subplot(5,1,4)
plot_fft(y_Poly0,128,1,1);title('直接滤波后频谱')
subplot(5,1,5)
plot(real(y_Poly0));title('直接滤波后数据')
 
 
figure(2)
subplot(4,1,1)
plot(y_sample);title('采样数据');
subplot(4,1,2)
plot_fft(y_sample,512,1,10);title('采样数据频谱');axis([-inf,inf,-100,5]);
subplot(4,1,3)
plot(real(mix_out));title('混频数据（I路）');
subplot(4,1,4)
plot_fft(mix_out,512,1,10);title('混频数据频谱');axis([-inf,inf,-100,5]);
 

截断函数使用round；

饱和函数：

function y=sat(x,bw,s)
%this function implement saturation process
% x:input signal, bw:bitwidth of x, s:1-signed,0-unsigned
if s == 1
   x(x>(2^(bw-1)-1))=2^(bw-1)-1;
   x(x<(-2^(bw-1)))=-2^(bw-1);
else
   x(x>(2^bw-1))=2^bw-1;
   x(x<(-2^bw))=-2^bw;
end
y = x;
end

频谱绘图：

function []=plot_fft(y_in,fs,db_en,p)%p可以在低点数下提高分辨率
N=length(y_in);
y_fft=abs(fftshift(fft(y_in,N*p)));
if db_en
    y_fft=20*log10(y_fft./max(y_fft));
else
    y_fft=(y_fft./max(y_fft));
end
plot(-fs/2:fs/N/p:fs/2-fs/N/p,y_fft);
end

滤波+抽取函数：

function [y_out]=FIR_deci(y_in,coe_in,D,bw)
y_in_conv=conv(coe_in,y_in);            % 128* s[16,15]*s[14,0] => s[30+7,15]
y_in_rnd = round(y_in_conv./2^bw); %去除小数 s[22,0]
y_in_sat = sat(y_in_rnd,16,1);  %去除6位整数 s[16,0]
y_out=y_in_sat(1:D:end);
 
end

滤波器设计：

RTBW = 100e6; % Real-time bandwidth = 20KHz
Fs = 512e6; % ADC sample rate
Fs_Eq = RTBW * 2.48; % Equivalent sample rate of the DDC I/Q outputs
DDCDecimationRate = Fs/Fs_Eq; % DDC decimation rate
 
InputWordLengths = 16; % Input word length of DDC (ADC output)
InputFracLengths = 15; % Input fraction length of DDC (ADC output)
CoefWordLengths = 16; % The word length of filter coefficients 
CoefFracLengths = 15; % The fraction length of filter coefficients
%% Polyphase decimator 0 Design (Decimation Rate = 4)
Fs_PolyDecim0 = Fs;
D_PolyDecim0 = 4;
N_PolyDecim0 = 127;
Fp_PolyDecim0 = RTBW / Fs_PolyDecim0;
Ap_PolyDecim0 = 0.05;
Ast_PolyDecim0 = 130;
 
H_PolyDecim0 = fdesign.decimator(D_PolyDecim0,'lowpass','N,Fp,Ap,Ast',...
               N_PolyDecim0,Fp_PolyDecim0,Ap_PolyDecim0,Ast_PolyDecim0);
Hd_PolyDecim0 = design(H_PolyDecim0,'SystemObject',true);
 
% fvtool(Hd_PolyDecim2)
 
% Creates a quantizer object
q_PolyDecim0 = quantizer([CoefWordLengths CoefFracLengths],'Mode','fixed',...
              'RoundMode','round','OverflowMode','saturate'); 
% Converts numeric matrix x into a hexadecimal string returned in y. 
% The attributes of the number are specified by the quantizer object q.
 
Poly0_coe=Hd_PolyDecim0.Numerator;
%生成系数文件
for k=1:4
PolyDecim0_Coef = num2hex(q_PolyDecim0,Poly0_coe(k:4:end));
[NumPolyDecim0Coef,~] = size(PolyDecim0_Coef);
 
PolyDecim0CoefFileName = strcat('Poly_phase_',num2str(k,'%d'),'.coe');
 
fileId = fopen(PolyDecim0CoefFileName,'w');
fprintf(fileId,'Radix = 16;\n');
fprintf(fileId,'Coefficient_Width = 16;\n');
fprintf(fileId,'CoefData = %s,\n', PolyDecim0_Coef(1,:));
for i = 2 : NumPolyDecim0Coef - 1
    fprintf(fileId,'%s,\n',PolyDecim0_Coef(i,:));
end
fprintf(fileId,'%s;\n', PolyDecim0_Coef(NumPolyDecim0Coef,:));
fclose(fileId);
end

5/16 

滤波器系数求和为1；绝对值求和在1-2；因此输入的s[14,0]与系数加权求和并去除15位小数后，得到的输出范围依旧为14位。

因为考虑到系数绝对值求和大于1，输出应该扩一位，而数据分为了IQ两路，两路分别有一半多项数据为0，输出位宽应该减一位，因此在不保留小数情况下，最后得到输出为两路s[14,0]数据。要输出16位数据，应该保留两位小数。

修改后的结构如下：