全国大学生电子设计竞赛2017年E题 自适应滤波器

本题拟采用实时采样的方式来进行自适应滤波。

学过数字信号处理都知道，频率分辨率是采样率除以采样点数，而这道题提高部分要求10Hz的分辨率，经过前期的分析，我们采用采样率为4MHz，那么就是说我们需要采样到至少400k个点才能进行一次自适应，而采样这么多点需要的时间是1/10 = 0.1s，所以每进行一次自适应就需要花费0.1s的时间。
而提高部分又要求要在1s内完成，所以我们大概可以进行7到8次的滤波。

经过matlab仿真不断移相噪声相减得到的相减波形能量如下结果

噪声是正弦波时，移相的能量变化图

噪声是三角波时，移相的能量变化图

噪声是方波时，移相的能量变化图

通过前期的matlab能量仿真我们发现，能量都可以近似看成周期性变化的单调函数。为了分析方便我这里直接把几种不同的噪声移相产生的能量都看成三角波。

算法思路

先找到一个这样的A点，该点能量斜率为负（即能量是随着移相点数增加而减小的），然后再在这个点之后找到能量上升的一个点，然后找到AB横坐标中点C，选取AB中能量小的那个和C组成新的能量对，继续找中点，直到迭代到最低点结束迭代。

下面是matlab算法验证

fa=10010;        %设置波形频率 
fb=10000;
Fs=16000000;        %设置采样频率
L=1600000 + 1600;         %数据长度
%=============产生输入信号==============%
t=0:1/Fs:(1/Fs)*(L-1);


A = sin(2*pi*fa*t+pi/3);
B = square(2*pi*fb*t);


D = sin(2*pi*fa*t+pi/3+pi)+square(2*pi*fb*t+pi);
%


%     for i = 1:1600
%     E = D(i:L-1600+i-1) - B(1:L-1600);
%     energy(i) = sum(E.^2);
%         
%         
%     end

%



%由于能量曲线的频率是根据噪声的周期来的，噪声的最大频率是100k，我们取两个大点的间隔的时候不能超过噪声的半个周期，所以最大的间隔点数小于40

%先找到第一个斜率下降的点
x_1 = 1;
x_2 = x_1 + 1;
x_3 = x_1 + 41;
x_4 = x_3 - 1;
energy_x_1 = cal_energy(D,B,x_1,L);
energy_x_2 = cal_energy(D,B,x_2,L);
energy_x_3 = cal_energy(D,B,x_3,L);
energy_x_4 = cal_energy(D,B,x_4,L);
while(slope(energy_x_2,energy_x_1,x_2,x_1) >= 0 || slope(energy_x_3,energy_x_4,x_3,x_4) <= 0)
    x_1 = x_1 + 41;
    x_2 = x_1 + 1;
    x_3 = x_1 + 40;
    x_4 = x_3 - 1;
    energy_x_1 = cal_energy(D,B,x_1,L);
    energy_x_2 = cal_energy(D,B,x_2,L);
    energy_x_3 = cal_energy(D,B,x_3,L);
    energy_x_4 = cal_energy(D,B,x_4,L);
end

     x_1 = x_2;
    x_mid = x_2 + 20;
    x_2 = x_2 + 40;
    energy_x_1 = cal_energy(D,B,x_1,L);
    energy_x_2 = cal_energy(D,B,x_2,L);
    
    if(energy_x_1 < energy_x_2)
       % x_1 = x_1;
        x_2 = x_mid;
    elseif(energy_x_1 > energy_x_2)
        x_1 = x_mid;
       % x_2 = x_2;
    else
        output = (x_1 + x_2)/2;
    end
    x_mid = (x_1 + x_2)/2;

%找到那个第一个斜率下降的点之后，开始进行二分搜索
for i = 1: 5
    energy_x_1 = cal_energy(D,B,x_1,L);
    energy_x_2 = cal_energy(D,B,x_2,L);
    if(energy_x_1 < energy_x_2)
       % x_1 = x_1;
        x_2 = x_mid;
    elseif(energy_x_1 > energy_x_2)
        x_1 = x_mid;
       % x_2 = x_2;
    else
        output = (x_1 + x_2)/2;
    end
    x_mid = (x_1 + x_2)/2;
end
   output = uint32((x_1 + x_2)/2 + 0.5);
    
E_lvbo = D - square(2*pi*fb*(t-double(output)./Fs));
figure;
plot(t,E_lvbo);  
Hd = LPF;
fil_E = filter(Hd,E_lvbo);
figure;
plot(t,fil_E);
    
    

其中slope函数如下

function y = slope(y2,y1,x2,x1)
    y = (y2 - y1)./(x2 - x1);

end

计算能量函数cal_energy如下

function y = cal_energy(D,B,i,L)
    i = uint32(i);
    E = D(i:L+i-1 - 1600) - B(1:L-1600);
    y = sum(E.^2);  
end

试题实现

第一部分：加法器

加法器采用硬件加法器，这个一个运放几个电阻就可以搭好，这里就不再赘述。

第二部分：移相器

移相器，由于硬件移相器对方波等存在高次谐波的信号移相效果略差，所以采用数字移相的方式来实现这部分内容。

这里我移相的方法采用硬件移相的方法，就是利用FPGA内部的FIFO存储数据以延时。所以我如果有两个信号，我对其中一个信号的延时即是两个信号的相位差发生了改变。

根据这个原理我们知道，FIFO存储的话一定是延时，所以在后续的自适应滤波器里面，我们要延时的信号也是B信号。

第三部分：自适应滤波器

MCU与FPGA部分实现

MCU要实现的函数：
1、发送延时量
2、计算斜率
3、判断最佳能量点

FPGA要实现的module
1、SPI
2、Sampling
FPGA要干的就是延时，然后采样计算能量发送，其他的没什么活

采样率8M，要分辨10Hz的信号，那么就是至少要采样8M/10个点，即总共要计算800k次，而每采样一次将花费大约0.1s的时间。

当采样完毕后将计算好的能量发送给MCU进行下一次的计算。