博客

clc;
close all;

Fs = 500;                                                     %设置采样频率
t = 0:1/Fs:3;  
t = t';
Size_t = size(t,1);
F1 = 7;
F2 = 13;
F3 = 23;
F4 = 50;
SIR = -100;            %信干比 Unit:dB

Signal = 10^(SIR/20)*(sin(2*pi*F1*t) + 0.5*sin(2*pi*F2*t) + 0.25*sin(2*pi*F3*t)); %生成信号
noise = 0.95*sin(2*pi*F4*t + pi/4);
Signal_noise = Signal + noise;                               %加入50Hz工频干扰
     
%RLS最小二乘法  仿真效果好于LMS

M = 2;                         %定义FIR滤波器阶数
lamda = 0.90;                  %定义遗忘因子     
Signal_Len = Size_t;           %定义信号数据的个数
I = eye(M);                    %生成对应的单位矩阵
c = 0.01;                      %小正数 保证矩阵P非奇异
y_out = zeros(Signal_Len,1);   %滤波器输出
Eta_out = zeros(Signal_Len,1); %误差输出
w_out = zeros(Signal_Len,M);   %系数输出
for i=1:Signal_Len
    %输入数据
    if i == 1                 %如果是第一次进入
        P_last = I/c;         %初始化滤波器系数
        w_last = zeros(M,1);  %初始化信号向量
    end
    d = Signal_noise(i);                         %输入新的期望信号
    x = [sin(2*pi*F4*(i-1)/Fs)
         cos(2*pi*F4*(i-1)/Fs)];                 %输入新的信号矢量
    %算法正体
    K = (P_last * x)/(lamda + x'* P_last * x);   %计算增益矢量
    y = x'* w_last;                          %计算FIR滤波器输出
    Eta = d - y;                             %计算估计的误差
    w = w_last + K * Eta;                    %计算滤波器系数矢量
    P = (I - K * x')* P_last/lamda;          %计算误差相关矩阵
    %变量更替
    P_last = P;
    w_last = w;
    %滤波结果存储
    y_out(i) = y;
    Eta_out(i) = Eta;
    w_out(i,:) = w';
end

figure;
subplot(3,1,1);
plot(t,Signal);
title('原始信号');
xlabel('时间t/s');
subplot(3,1,2);
plot(t,Signal_noise);
string = ['加入50Hz工频干扰的信号,信干比SIR = ',num2str(SIR),'dB'];
title(string);
xlabel('时间t/s');
subplot(3,1,3);
plot(t,Eta_out,'b',t,Signal,'r');
%开始不拟合 
title('输出误差');
xlabel('时间t/s');
axis([0,3,-2*10^-5,2*10^-5]);%