利用自相关算法去除50Hz噪声（构想）

软件50hz陷波，脱离传统陷波策略，考虑下自相关滤波：提取出50hz成分，然后减去这个成分。

把空信号（只有白噪声和50hz）的信号进行自相关，得到自相关值和延迟时间，从而得到50hz的频率，在ifft 还原50hz信号。最后在有数据的信号中减去这个50hz信号
自相关，也叫序列相关，是一个信号与其自身在不同时间点的互相关。非正式地来说，自相关是对同一信号在不同时间的两次观察，通过对比来评判两者的相似程度。自相关函数就是信号x(t)和它的时移信号x(t-τ)的乘积平均值。它是时移变量τ的函数。

在这里附上自相关的公式 😄

step1：在matlab中建立各种频率信号，如1Hz的低频信号、100Hz的高频信号、50Hz工频干扰以及白噪声（其中1Hz和100Hz被视为信息信号）

clc;
clear;
close all;

n = 3000;
fs = 1000;
t = (0:n-1)/fs;
ff = fs*(0: n/2)/n;

f0 = 1;
f1 = 100;
fn = 50;
N = 2000;%横轴显示范围
T = t(1:N);

x0 = 5*sin(2*pi*f0*t);
x1 = 3*sin(2*pi*f1*t);
xn = sin(2*pi*fn*t);

white_noise = 0.1*randn(size(t)); %白噪声

y_null = xn + white_noise;%表示没有或直流信号输入
y_low = xn + x0 + white_noise;%表示在输入信号为低频信号
y_high = xn + x1 + white_noise;%表示输入信号为高频信号
y_high_clean = x1 + xn;
y_mix = xn + x1 + x0 + white_noise;%表示输入信号为混合频率信号

step2:对空白信息信号（仅仅包括白噪声和50Hz工频干扰）进行自相关运算

[R, tau] = xcorr(y_null,100,'coeff');

plot(tau,R);
xlabel('滞后');
ylabel('自相关函数');

不难看出，即使我改变白噪声的权重，在延时时间tau=0之后，依然表现出明显的周期信号，这个就是作为噪声的周期干扰（50Hz工频干扰）

step3:由于得到的图像是离散的所以为了关注整体波动趋势而过滤小幅度高频信号的波动。我想到了包络线来关注平台的变化。（三次样条插值）

下面的lift应该是right才对，sry！！！👀

R_lift = R(101:end); %因为自相关函数是对称的，因此只取右边. 
R_lift = R_lift';

c=zeros(size(R_lift));
d=zeros(size(R_lift));
s=zeros(size(R_lift));
decrease=zeros(size(R_lift)); 
increase=zeros(size(R_lift)); 
spc=100;           % spc=sample/cycle
n1=0:360/spc:360;
n1=n1';

dp=diff(R_lift,1)/(360/spc); 
for i=1:size(dp,1)-1
 if dp(i)*dp(i+1)<0
  c(i)=n1(i+1); %寻找驻点横坐标
  d(i)=R_lift(i+1);  %寻找驻点纵坐标
 end
end 
c(find(c==0))=[]; 
d(find(d==0))=[];

figure(2);
plot(n1,R_lift); 
hold on;

%%
for i=1:size(c,1)
 s(i)=find(n1==c(i)); %驻点横坐标的序列号 
end
s(find(s==0))=[]; % 序列更紧凑，下同
j=1;
for i=1:size(c,1)-1
 if R_lift(s(i+1))>R_lift(s(i)) %寻找下降区间（s,decrease为序列号） 
     decrease(j)=s(i); %i~i-1为下降区间（s,decrease为序列号）
        increase(j)=s(i+1); %i+1~i为上升区间（s,increase为序列号）
  j=j+1; %i~i+1为上升区间（序列号） 
 end
end
decrease(find(decrease==0))=[]; 
increase(find(increase==0))=[]; 

peak_n1=n1(increase);
peak_R_lift=R_lift(increase)
% 由于寻找极大值会漏掉原始数据的第一个和最后一个点，从而造成包络线在初末位置不完整，所以需要进行额外处理。此处由于p1是周期序列，所以对原始序列p1最后一个值赋值到peak_p1序列的第一个值和最后一个值即可。
peak_R_lift(2:size(peak_R_lift,1)+1)=peak_R_lift(1:end)
peak_R_lift(1)=R_lift(end)
peak_R_lift(size(peak_R_lift,1)+1)=R_lift(end)
% n1 同理
peak_n1(2:size(peak_n1,1)+1)=peak_n1(1:end)
peak_n1(1)=0
peak_n1(size(peak_n1,1)+1)=n1(end)
scatter(peak_n1,peak_R_lift) % 绘制极大值点
%对新的数据进行插值，此处三次样条将极大值序列插值回原始数据长度
peak_pp=interp1(peak_n1,peak_R_lift,n1,'spline');
peak_R_lift=peak_pp;
figure(2);
plot(n1,peak_R_lift);

图中黄色的线已经画出包络线，并且标记出了极大值点，好了，下一步就是算出每个极大值的横坐标之间的差值了。

step4:我采用了计算标记出所有的极大值点来求算术平均来得到周期（取整了）

%% 计算极大值的平均
D_val = 0;
D_sum = 0;
time = 0;
for i = 2: size(increase)
    D_val = increase(i) - increase(i-1);
    D_sum = D_sum + D_val;
    time = time + 1;
end

period = ceil(D_sum/time)/1000;
period_freq = 1/period;

margin = 1;  %这个margin是为下面在频域中截取做预留值

step5:下面就是在频域里去除掉我们不要的部分，得到我们需要的部分，并且还原特定频率的信号成分

df=fs/n;            % 频率间隔
f=(0:1:n/2)*df;     % 频率刻度
f=f';               % 为了方便查看, 将行向量 f 转置成列向量

figure(1)

subplot(3,2,[1,2]);
plot(t, y_null);                    

subplot(3,2,[3, 4])
plot(ff, Y)          
xlabel('频率/Hz')
ylabel('幅值/bar')

X(1:period_freq*(n/fs)+1-margin)=0; %去除除了50频以外所有频率
X(period_freq*(n/fs)+1 + margin: n - period_freq*(n/fs)+1-margin)=0;  %去除除了50频以外所有频率
X(n - period_freq*(n/fs)+1+margin: end)=0;
IX = ifft(X);

subplot(3,2,[5,6])
res = IX+mean(y_null);
plot(T, res(1:N));

%  弥补直流分量
ff_50 = IX + mean(y_null);

太好了，我们提取到了干净的50Hz成分。✌️

step6:接下来我们试试在有信息信号中减去这个50Hz工频干扰的效果怎么样

figure(3);
hold on;
plot(T, y_low(1:N), 'r');
plot(T, y_low(1:N) - ff_50(1:N), 'b');

summary: 这个方法看起来很有意思，但是，其实存在很多的问题：比如，50Hz在做减法的时候，没有考虑相位的问题，当空信息信号中的50Hz相位与载信息信号中的50Hz相位相差180°时不但不会起到正向的效果，反而，会使50Hz的幅值翻倍😣

more: 我想如果这个方法被用在上位机中，集成上可以调节50Hz相位的功能，通过观察SNR或许可以起到不错的效果。