当前位置:   article > 正文

数字信号处理学习:基于CIC滤波器的正交检波_cic数字

cic数字

基于CIC滤波器的正交检波

 

Abstract:

本文主要介绍一种基于CIC滤波器的正交检波结构,前部分主要介绍CIC滤波器的实现,并结合MATLAB进行图文分析;后半部分简单介绍正交检波的一些个人理解。

Keywords:

CIC滤波器、抽取、混叠、FFT、欠采样

 

CIC滤波器初步介绍

Xilinx的FPGA的DSP IP中有一个模块就是CIC滤波器,这个模块可以不用乘法器资源(硬件资源消耗少)来实现低通滤波、内插、抽取:在作为低通滤波器使用时,因为没有用乘法器,所以FPGA内部的delay就可以降到最低,以实现更大的数据吞吐率。同时基于CIC滤波器可以较好地实现数据抽取,已达到降采样率的目的(相同的滤波器特性情况下,降低采样率可以进一步的减少FIR滤波器的阶数,也就能够进一步的减少FPGA的资源使用)。结合这两个特点,CIC滤波器很适合放置在高采样率的ADC的后面,1实现初步滤波,2利用抽取降低采样率,使得后级更加陡峭的滤波器更加容易应用。下面是CIC滤波器的Z域传函:

Figure1

下面是CIC的实现方式:

Figure2

左边是积分部分,共N级,时钟域为降采前(假设R为2)的时钟,如果前级连接一个10MHz的ADC,那么Integrator部分的计算时钟就是10MHz。右边部分是Comb部分,共N级,每一级有M个延时单元,时钟域是降采后的时钟,接上面的例子,就是5MHz。

CIC滤波器仿真

下面给出基于MATLAB的源代码:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function yn = CIC_Simulat(R,M,N,xn )

yn = zeros(1,fix(length(xn)/R));

IntegratorBuffer = zeros(1,N);

IntegratorIn = 0;

IntegratorOut = 0;

CombBuffer = zeros(N,M);

CombIn = 0;

CombOut = 0;

for i_Sample = 1 : length(xn)

    IntegratorIn =xn(i_Sample);

    %without delay, get IntegratorOutput at first

    %then update the IntegratorBuffer

    fori_Integrator = 1 : N

       IntegratorOut = IntegratorIn + IntegratorBuffer(i_Integrator);

       IntegratorBuffer(i_Integrator) = IntegratorOut;

       IntegratorIn = IntegratorOut;

    end

    %case downsampling, input data to comb from every R IntegratorOutput

    ifrem(i_Sample,R) == 0

        CombIn=  IntegratorOut; % CombIn is Nth Comb's input sigal

        fori_Comb = 1 : N

            CombOut= CombIn - CombBuffer(i_Comb,end); % CombOut is Nth Comb'soutput sigal

            %thenupdate the Nth Comb's CombBuffer

           CombBuffer(i_Comb,2:end) = CombBuffer(i_Comb,1:end-1);

           CombBuffer(i_Comb,1) = CombIn;

            CombIn= CombOut;% reset CombIn as next N+1th Comb's input

        end %end for

        %at last we get output

       yn(fix(i_Sample/R)) = R * CombOut/((R*M)^(N));     

    end%end if

end%end for

 

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

基于上面的脚本,来简单的做个实验:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 该脚本主要用来展示,LNFFT的影响

% L越长,也就是实际用来做FFT的采样数据点数越多,那么加窗效应就会减弱,主瓣的宽度就会越小,功率泄露就会减弱(这里特指矩形窗)

% N越长,就是在实际L的基础上,添加的0采样点的个数越多,就会提高频谱的分辨率,也就是最小能够分辨的频率

fftN = 1024; % fft number

sF = 100; %sample frequence

sL = 256; %sample Length

 

 

xF1 = 1; %x(n) frequence 1

xF2 = 37; %x(n) frequence 2

xA1 = 1;

xA2 = 1;

xn = xA1 * sin(2*xF1*pi*(0:(1/sF):((sL-1)/sF)))+...

     xA2 *sin(2*xF2*pi*(0:(1/sF):((sL-1)/sF)));%Sample freq = sF, length = sL

xF = fft(xn,fftN);

figure

subplot(2,4,1);

plot((((-fftN/2):(fftN/2-1))/fftN)*sF,abs(fftshift(xF))); 

xlabel('Normalized frequency') ;

ylabel('Amplitude') ;

title('Input Signal(FFT)');

subplot(2,4,5);

plot((((-fftN/2):(fftN/2-1))/fftN)*sF,abs(fftshift(xF)));

xlabel('Normalized frequency') ;

ylabel('Amplitude') ;

title('Input Signal(FFT)');

 

 

%-----------------------Average-----------------------------%

h1L = 10; %h(n) Length

h1n = ones(1,h1L); %get h1(n)

subplot(2,4,2);

plot(h1n,'Ro');

xlabel('Time') ;

ylabel('Amplitude') ;

title('Impulse Response of Average');

 

 

h1F = fft(h1n,fftN)/h1L;

subplot(2,4,3);

plot((((-fftN/2):(fftN/2-1))/fftN)*sF,abs(fftshift(h1F)));

xlabel('Normalized frequency') ;

ylabel('Amplitude') ;

title('Magnitude Response Of Moving Average N=10 (FFT)');

 

y1n = conv(xn,h1n)/h1L;%get y1n

y1F = fft(y1n,fftN); 

subplot(2,4,4);

plot((((-fftN/2):(fftN/2-1))/fftN)*sF,abs(fftshift(y1F)));

xlabel('Normalized frequency') ;

ylabel('Amplitude') ;

title('Output Signal Of Moving Average Filter (FFT)');

%------------------------------------------------------------%

 

 

 

%-----------------------CIC-----------------------------%

% Implementing Cascaded Integrator Comb filter with the

% integrator section following the comb stage

% RMN来控制CIC滤波器的行为,其中

% 实际上RMN都会影响到带宽,最先确定的一般是R,就可以通过设置MN来实现CIC滤波器的幅度响应

% R来控制Downsampling,也即是下采的比例,一般固定不变(会影响衰减的幅度和旁瓣的宽度)

% M来控制带宽,(会影响衰减的幅度和旁瓣的宽度)

% N来控制CIC滤波器的阶数,主要用来控制衰减的的幅度

R = 1;

M = 10;

N = 1;

%理论CIC的冲激响应,冲激响应FFT,信号输出FFT

 

%理论CIC的冲激响应

z=tf('z');

h3Theory = ((1-z^(-M*R))/(1-z^(-1)))^N;

h3nTheory = impulse(h3Theory)/((R*M)^(N));

%h3nTheory = step(h3Theory);

%实际CIC的冲激响应

DeltaXn = [1 zeros(1,length(h3nTheory)-1)];

%DeltaXn = ones(1,length(h3nTheory));

h3nReal = CIC_Simulat(R,M,N,DeltaXn);

 

 

subplot(2,4,6);

hold on;

plot(h3nTheory,'G*');

plot([R:R:R*length(h3nReal)],h3nReal,'R.');

xlabel('Time') ;

ylabel('Amplitude') ;

title('Impulse Response of CIC');

hold off;

 

%理论CIC的冲激响应的FFT

h3TheoryF = fft(h3nTheory,fftN);

%实际CIC的冲激响应的FFT

h3RealF = fft(h3nReal,fftN);

subplot(2,4,7);

hold on;

plot((((-fftN/2):(fftN/2-1))/fftN)*sF,abs(fftshift(h3TheoryF)),'G-');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF/R),abs(fftshift(h3RealF)),'R-');

xlabel('Normalized frequency') ;

ylabel('Amplitude') ;

title(' Magnitude Response Of CIC (FFT)');

hold off;

 

 

y3n=CIC_Simulat(R,M,N,xn);

y3F = fft(y3n,fftN); 

subplot(2,4,8);

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF/R),abs(fftshift(y3F)));

xlabel('Normalized frequency') ;

ylabel('Amplitude') ;

title('Output Signal Of CIC Filter (FFT)');

%------------------------------------------------------------%

 

上面的代码主要是做了一个对比实验,因为当R=1,N=1,M=10时,CIC滤波器就是实现了移动平均滤波器。将CIC滤波器以及移动平均滤波器的冲激响应和冲激响应的FFT进行对比,结果如下图所示:

Figure 3 移动平均滤波器:输入信号的频谱、冲激响应、冲激响应幅频曲线、滤波后的输出信号的频谱

Figure 4 CIC滤波器(R=1、M=10、N=1):输入信号的频谱、冲激响应、冲激响应幅频曲线、滤波后的输出信号的频谱

经过对比,两者确实具备相同的冲激响应和幅频特性,结果输出信号自然是一样。

现将R=2,M=10,N=2,,得到下图:

Figure 5 CIC滤波器(R=2、M=10、N=2):输入信号的频谱、冲激响应、冲激响应幅频曲线、滤波后的输出信号的频谱

特别需要注意的是上图左2,红色点和绿色点没有重合,同时绿色点数量也是红色点的2倍。这是因为R=2,这里有一个降采样率(抽取)的过程,那自然冲激响应的有效点数也是原来的一半,幅度却是两倍(红线是CIC抽取的响应,绿点是没有降采的响应)。再看右2图,因为抽取的原因,红线的频率范围只有绿线的1半,但是频率重合部分的幅频响应是一致的。所以说同样的传函形式:

 

 

 

可以用降采(抽取)的算法实现(红线),也可以不降采(绿线)来实现,最终幅频响应是一样的。

再回过头对比figure1和figure3,可以明显的看到,高阶CIC高频部分衰减还是可以接受的。但是CIC抽取滤波器存在两个问题:

1.     抽取时会发生频谱折叠(或说频谱复制),有可能会发生混叠;

做如下修改:

xF2 = 42.5; %x(n) frequence 2

R =2;

M =10;

N =1;

 

此时两个滤波器的实际输出如下:

Figure7

可以看到,原来的频率应该是1Hz和42.5Hz,但是CIC滤波器后变成了1Hz和7.5Hz,这是因为在降采后发生了频率复制(折叠),那原来在42.5Hz的频率向左移动50Hz(降采后的采样率)-7.5Hz,同理-42.5Hz移动到了7.5Hz。这就是混频,一般发生在欠采样(抽取后采样率不够大)的情况下!

2.     带通部分平整度不够;

Figure8 CIC通带部分展开后,却是不够平整

前者可以通过适当选取采样频率进行规避,后者可以在CIC滤波器后面接一个补偿器,将通带部分的幅频响应变得更加平整,如下图所示:

Figure9 经过补偿的CIC滤波器,见Xilinx的技术手册-pg140

小结:

上面主要介绍了CIC滤波器的特点,接下来介绍下应用。

正交检波介绍:

Figure10 图片来自《用于低场磁共振成像的直接射频采样数字接收机》

具体的理论可以看信号与系统等书籍,这里先给出MATLAB的仿真代码:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

R = 2;

M = 10;

N = 1;

 

sN = 5120;             %sample number per fft

fftN = 10240;        %fft number

sF1M = 1e6;            %sample qrequence

sF5M = 5e6;            %sample qrequence

fc = 1.11e6;         

                       %R=2,N=1,M=10测试环境

                       %混频1 载波中频 1.4625e6(fc-fs<0.5fs) :左上图两根蓝线就是用5M去采样获得的FFT,红线就是用1M去采样获得的FFT

                       %右边的红线0.4625就是右边的蓝线1.4625频移1Mfs)得到,左边的红线-0.4625是由左边的蓝线-1.4625频移1M得到

                       %实际上在左边的红线的左边,还有一根很被隐藏了的红线-0.5375,以及右边有一根红线0.5375。因为不落在±0.5之内

                       %但是在用NCO进行频谱搬移后,会发现这根线-0.5375的线移动到±0.075处,也就是移动到了实际有效带宽范围,如左下图所示

                       %这个时候看右下图的检波信号,会发现有大量的高频谐波在上面,正是这个原因!1.5>fc>1.375会有这个问题

                       

                       %混频2

                       %载波中频为1.251.375之间时,虽然上述干扰谱线在经过NCO移频后已经在±0.5之外了,但是如果R大于等于2,就会有进一步频谱幅值

                       %但是一般这时的干扰谱线已经经过充分的滤波,所以即使混叠了,影响也不大!

                       

                       %实际上都需要结合中频带宽,中频中心频率,采样频率,降采样率共同来看是否有混叠

 

 

fa = [0.0022e6 0.0032e6 0.0042e6 0.0052e6];

fp = 0.0013e6;

fnco =rem(fc,sF1M);

 

t1M=0:1/sF1M:(sN-1)/sF1M;

t5M=0:1/sF5M:(sN-1)/sF5M;

 

 

an1M =cos(2*pi*fa(1)*t1M)+sin(2*pi*fa(2)*t1M)+cos(2*pi*fa(3)*t1M)+sin(2*pi*fa(4)*t1M);

an5M =cos(2*pi*fa(1)*t5M)+sin(2*pi*fa(2)*t5M)+cos(2*pi*fa(3)*t5M)+sin(2*pi*fa(4)*t5M);

 

 

 

pn1M = pi*cos(2*pi*fp*t1M);

pn5M = pi*cos(2*pi*fp*t5M);

 

xn1M = an1M .* sin(2*pi*fc*t1M + pn1M);

xn5M = an5M .* sin(2*pi*fc*t5M + pn5M);

x1MF = fft(xn1M,fftN)/sN;

x5MF = fft(xn5M,fftN)/sN;

 

ncoI1M = cos(2*pi*fnco*t1M);

ncoQ1M = sin(2*pi*fnco*t1M);

ncoI1MF = fft(ncoI1M,fftN)/sN;

ncoQ1MF = fft(ncoQ1M,fftN)/sN;

 

 

sIn1M = xn1M .* ncoI1M;

sQn1M = xn1M .* ncoQ1M;

sIn1MF = fft(sIn1M,fftN)/sN;

sQn1MF = fft(sQn1M,fftN)/sN;

 

sIcic1RM = CIC_Simulat(R,M,N,sIn1M);

sQcic1RM = CIC_Simulat(R,M,N,sQn1M);

sIcic1RMF = fft(sIcic1RM,fftN)/sN;

sQcic1RMF = fft(sQcic1RM,fftN)/sN;

 

 

sn = sIcic1RM + i * sQcic1RM;

sMag = abs(sn);

 

 

%理论CIC的冲激响应

z=tf('z');

h3Theory = ((1-z^(-M*R))/(1-z^(-1)))^N;

h3nTheory = impulse(h3Theory)/((R*M)^(N));

%实际CIC的冲激响应

DeltaXn = [1 zeros(1,length(h3nTheory)-1)];

%DeltaXn = ones(1,length(h3nTheory));

hnReal = CIC_Simulat(R,M,N,DeltaXn);

hRealF = fft(hnReal,fftN);

subplot(2,2,1);

 

 

 

hold on;

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M),abs(fftshift(x1MF)),'R-');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF5M),abs(fftshift(x5MF)),'B--');

xlabel('Frequence') ;

ylabel('Amplitude') ;

title('Real(Blue) vs Sample(Read) Magnitude');

hold off;

 

subplot(2,2,2)

hold on;

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M),abs(fftshift(x1MF)),'R--');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M),abs(fftshift(ncoI1MF)),'G--');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M),abs(fftshift(sIn1MF)),'B--');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M/R),abs(fftshift(sIcic1RMF)),'k-');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M/R),abs(fftshift(hRealF)),'y-');

xlabel('Frequence') ;

ylabel('Amplitude') ;

title('Sample(Red) vs NCOI(Green) vs SI(Blue) vsSICIC(Black)Magnitude');

hold off;

 

subplot(2,2,3)

hold on;

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M),abs(fftshift(x1MF)),'R--');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M), abs(fftshift(ncoQ1MF)),'G--');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M),abs(fftshift(sQn1MF)),'B--');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M/R),abs(fftshift(sQcic1RMF)),'k-');

plot((((-fftN/2):(fftN/2-1))/fftN)*(sF1M/R),abs(fftshift(hRealF)),'y-');

xlabel('Frequence') ;

ylabel('Amplitude') ;

title('Sample(Red) vs NCOQ(Green) vs SQ(Blue) vsSQCIC(Black) Magnitude');

hold off;

 

subplot(2,2,4)

hold on;

plot(((0:(sN-1))/sF1M),xn1M,'R-');

plot(((0:(sN-1))/sF5M),xn5M,'G-');

plot(((0:((sN/R)-1))/(sF1M/R)),sMag,'k--');

xlabel('Sample Time') ;

ylabel('Amplitude') ;

title('Real(Green) vs Sample(Red) vs |S|(Black)Magnitude');

hold off;

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

将载波(中频)fc设置为1.11e6,这个时候我们用两种采样率去采集这个波形,得到的图案是:

Figure11 红线是1M采样率采样后所做的FFT,蓝线是5M采样率的FFT

可以看到,用5MHz进行采样时能够获得正确的频谱,而用1MHz进行采样时,获得的就是经过复制后的频谱,其中原来蓝线1.11M处的跑到了0.11M(红线),而原来蓝线-1.11M的跑到了-0.11M。这是一个典型的欠采样,也就是用低采样率去采集中频信号,然后利用频谱的复制,得到可操作的频谱。

这时再用NCO生成一个余弦和正弦分别乘以ADC后的数据,可以得到如下SI和SQ频谱:

Figure12 SI:将ADC的数据(红)乘以余弦(绿)得到SI(蓝),同时经过CIC滤波器(幅频特性黄线)后变成SICIC(黑)

Figure13 将ADC的数据(红)乘以余弦(绿)得到SQ(蓝),同时经过CIC滤波器(幅频特性黄线)后变成SQCIC(黑)

通过上面两张图,可以看到ADC乘以NCO信号后,将红线分别挪到了蓝线处,其中一处蓝线在2.2M,另一处在0M(蓝线被黑线覆盖了)。然后经过CIC滤波器后进行滤波,将2.2M处的高次谐波滤掉,就得到了干净的SICIC或者SQCIC。

将SQCIC和SICIC的合并:

sn = sIcic1RM + i * sQcic1RM;

sMag = abs(sn);

最终得到了检波信号的时间域信息,如下图:

Figure 14 绿线是5M采样,红线是1M采样,黑色虚线是检波信号(代表包络,有延时)

Figure 15 放大后:绿线是5M采样,红线是1M采样,黑色虚线是检波信号(代表包络,有延时)

可以看到这种方式确实可以实现检波;

接下来分析一下特殊情况:

%R=2,N=1,M=10测试环境

混频1 载波中频 1.4625e6(fc-fs<0.5fs) :左上图两根蓝线就是用5M去采样获得的FFT,红线就是用1M去采样获得的FFT,右边的红线0.4625就是右边的蓝线1.4625频移1M(fs)得到,左边的红线-0.4625是由左边的蓝线-1.4625频移1M得到实际上在左边的红线的左边,还有一根很被隐藏了的红线-0.5375,以及右边有一根红线0.5375。因为不落在±0.5之内但是在用NCO进行频谱搬移后,会发现这根线-0.5375的线移动到±0.075处,也就是移动到了实际有效带宽范围,如左下图所示。这个时候看右下图的检波信号,会发现有大量的高频谐波在上面,正是这个原因!1.5>fc>1.375会有这个问题。

Figure16 这个图看不出异常

Figure17 问题来了

Figure 18 放大了,发现有谐波混叠进来了,而且恰好落在了旁瓣里,没有滤干净

 

下图可以看到,检波信号里面确实混杂了很多高频谐波

Figure 19

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/很楠不爱3/article/detail/698667
推荐阅读
相关标签
  

闽ICP备14008679号