数字滤波器的MATLAB和FPGA实现-读书笔记 (3/3)

第8章 变换域滤波器的FPGA实现

原理：

变换域滤波：

最常用的是频域处理，即将时域信号通过快速傅里叶变换（Fast Fourier Transform，FFT）变换到频域进行滤波处理，再将频域信号通过快速傅里叶反变换（Inverse Fast Fourier Transform，IFFT）变换成时域信号，从而完成滤波处理。

DFT存在的问题：

栅栏效应和序列补零，频谱泄漏和混叠失真，频率分辨率与DFT参数的选择

DFT和FFT算法计算量：

DFT算法，乘法次数与N^2成正比

FFT算法，(N / 2)log2N次复乘运算

变换域滤波技术：

变换域技术通过合适的变换把扩频信号、热噪声以及干扰组成的混合信号映射到另一个域处理，一般将NBI映射为类似于冲激的函数，将有用信号映射为与干扰正交，具有平坦特性的波形。这样就能在较彻底去除干扰的同时对有用信号的损伤较小。

MATLAB实现：

FFT算法MATLAB仿真：

仿真FFT参数对采用FFT算法分析信号频谱的影响。产生频率分别为f 1 =2 Hz、f 2 =2.05 Hz的正弦波合成信号，采样频率f s =10 Hz。根据式（8-8），要实现分辨两个单频信号的目的，DFT的序列长度必须满足N> 400。分别仿真3种情况下的FFT变换：（1）取x (n )的128点数据，计算FFT；（2）将128点x (n )以补零的方式加长到512点，计算FFT；（3）取512点x (n )，计算FFT。

f1=2; f2=2.05;            %单正弦信号的频率
fs=10;              %采样频率
%对 128 点时域序列进行 FFT 分析
N=128;              %FFT 分析的点数
n=0:N-1;
xn1=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);  %产生 128 点时域信号序列
XK1=fft(xn1);            %进行傅里叶变换，并进行归一化处理
MXK1=abs(XK1(1:N/2));
%对补零后的 512 点时域序列进行 FFT 分析
M=512;
xn2=[xn1 zeros(1,M-N)];        %在时域信号序列后补零
XK2=fft(xn2);            %进行傅里叶变换，并进行归一化处理
MXK2=abs(XK2(1:M/2));
%对 512 点时域序列进行 FFT 分析
n=0:M-1;
xn3=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);  %产生 128 点时域信号序列
XK3=fft(xn3);            %进行傅里叶变换，并进行归一化处理
MXK3=abs(XK3(1:M/2));
%绘图
subplot(321);
x1=0:N-1;
plot(x1,xn1);xlabel('n','fontsize',8);title('128 点 x(n)','fontsize',8);
subplot(322);
k1=(0:N/2-1)*fs/N;
plot(k1,MXK1);xlabel('f(Hz)','fontsize',8);title('128 点 xn 的 FFT 变换','fontsize',8);
subplot(323);
x2=0:M-1;
plot(x2,xn2);xlabel('n','fontsize',8);title('512 点补零 x(n)','fontsize',8);
subplot(324);
k2=(0:M/2-1)*fs/M;
plot(k2,MXK2);xlabel('f(Hz)','fontsize',8);title('512 点补零 xn 的 FFT 变换','fontsize',8);
subplot(325);
plot(x2,xn3);xlabel('n','fontsize',8);title('512 点 x(n)','fontsize',8);
subplot(326);
plot(k2,MXK3);xlabel('f(Hz)','fontsize',8);title('512 点 xn 的 FFT 变换','fontsize',8);

频域滤波器MATLAB仿真：

假设在伪码码率为1 MHz的扩频信号中有一个频率为100 kHz的单音干扰，且信干比为0 dB，信噪比为10 dB。采用256点FFT/IFFT变换，窗函数为汉明窗，50%重叠，门限因子θ =20。需要仿真出滤波前后信号的频谱及时域信号情况，同时需要将滤波器的输入数据量化为12比特数据，存放在外部文本文件中，作为后续FPGA仿真测试的输入源。

Fpoint=256;%fft 变换的点数为 256
poly=[1 0 0 0 1 1 1 1];         %产生 PN 码的多项式
%根据多项式产生伪码，同时对伪码进行 4 倍采样
%E8_2_PNCodeGenerator 为自定义的函数，请在光盘中查看函数代码
pn1=E8_2_PNCodeGenerator(poly);
window1=hamming(Fpoint)';      %产生 Hamming 窗
%生成相同长度的伪码 PN 及窗函数数据 Window
pn=pn1;
window=window1;
for i=1:9
        pn=[pn,pn1];
end
PN=pn(1:floor(length(pn)/Fpoint)*Fpoint);
N=length(PN);
for i=1:N/Fpoint-1
        window=[window,window1];
end
%产生功率为 1   W，频率为 100   kHz 的单音干扰信号
f=100000;              %信号频率为 100 kHz
Fs=4*10^6;            %采样频率为 4 MHz
t=0:1/Fs:(N-1)/Fs;
jam=sqrt(2)*sin(2*pi*f*t);
S=PN+jam;            %产生 SJR=0 dB 的测试信号
S2=[S(Fpoint/2+1:N),S(1:Fpoint/2)];    %对输入信号进行 50%重叠
WS=S.*window;           %加窗处理
WS2=S2.*window;          %加窗处理
%对两路加窗信号进行变换域滤波处理
Fout=zeros(1,N);
for n=1:N/Fpoint
        FWS=fft(WS((n-1)*Fpoint+1:n*Fpoint));
        FWS2=fft(WS2((n-1)*Fpoint+1:n*Fpoint));
        AWS=abs(FWS).*abs(FWS);
        AWS2=abs(FWS2).*abs(FWS2);
        MAWS=sum(AWS)/Fpoint
        MAWS2=sum(AWS2)/Fpoint;
        for i=1:Fpoint
                if AWS(i)>(20*MAWS)
                        FWS(i)=0;
                end;
                if AWS2(i)>(20*MAWS2)
                        FWS2(i)=0;
                end;
        end;
        iS(((n-1)*Fpoint+1):n*Fpoint)=ifft(FWS);
        iS2(((n-1)*Fpoint+1):n*Fpoint)=ifft(FWS2);
end
%对滤波后的两路数据进行重叠处理后输出
Fout=iS+[iS2(N-Fpoint/2+1:N),iS2(1:N-Fpoint/2)];
%绘图
Fs=fft(S,256);as=abs(Fs);
Fsout=fft(Fout,256);asout=abs(Fsout);
figure(1);              %绘制滤波前后信号的频谱图
subplot(211);plot(abs(as));legend('叠加干扰的原始数据频谱(sjr=0dB)');
subplot(212);plot(asout);legend('滤波后的原始数据频谱');
figure(2);              %绘制滤波前后信号的时域波形图
p=1:200;
subplot(311);plot(p,PN(1:200)>0,'-',p,S(1:200)>0,'--',p,Fout(1:200)>0,'+');
axis([0 200 -0.2 1.2]);title('判决整形后的数据');xlabel('point');
legend('原始数据','滤波前数据','滤波后数据');
subplot(312);plot(S(1:200));legend('叠加干扰的原始数据');xlabel('point');
subplot(313);plot(Fout(1:200));legend('滤波后的原始数据');xlabel('point');
%将输入信号进行量化，并写入外部文本文件中，供 FPGA 测试仿真使用

FPGA实现：

 直接的变换域滤波器结构：

 基于模块复用的变换域滤波器结构：

第9章 解调系统滤波器的FPGA实现

数字接收机：

Costas环：