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]);%
基于RLS自适应滤波器
于 2022-03-02 18:09:29 首次发布