【离散数据处理】Fourier平滑法（零相移）

最新推荐文章于 2024-06-12 10:00:23 发布

sanfan66

最新推荐文章于 2024-06-12 10:00:23 发布

阅读量1k

点赞数 14

文章标签：算法信号处理傅里叶分析

本文链接：https://blog.csdn.net/sanfan66/article/details/136334410

版权

【离散数据处理】快速傅里叶变换与滤波。

针对离散信号经普通数字滤波器如巴特沃斯滤波器后产生相移问题，介绍零相移的Fourier平滑法，方法基于Fourier分析，属于频域全局平滑方法，由于Fourier平滑法不仅需要之前和当前数据，还需要往后的数据，仅适用于离线数据分析。

滤波方法比较

$\begin{aligned} y_1&=2 \sdot sin(2\pi f_1\sdot t) \nonumber\\ y_2&= sin(2\pi f_2\sdot t-90°) \nonumber\\ y&=y_1+y_2 \end{aligned}$
$f_1=1Hz、 f_2=5Hz$ 、 $f_s=100Hz$ ，想要将 $y$ 中的 $y_2$ 信号去除，Fourier平滑法和巴特沃斯低通滤波器比较如下：

滤波效果比较

可见Fourier平滑法能够有效去除较高频信号，并且相移很小。

详细代码

main.m

%
% Auther:  sanfan66
%
clc
clear all
close all

fs=100;
t=0:1/fs:5;
f1=1;
f2=5;
y1=2*sin(2*pi*f1*t);
y2=sin(2*pi*f2*t-pi/2);
y.y=y1+y2;

figure
plot(t,y1,t,y2,t,y.y)
legend('y1','y2','y')
xlabel('t')
ylabel('y')
xlim([t(1) t(end)]);
grid on;grid minor;
saveas(gcf, '1 y', 'png')

%% Fourier平滑法
fmax=3;
y.Fourier = FourierFil(fmax,fs,y.y);
%% 巴特沃斯滤波器
wp=2;%通带截止频率
ws=4;%阻带截止频率
rp=3;%通带最大衰减
rs=30;%阻带最小衰减

[n,wn]=buttord(wp/(fs/2), ws/(fs/2), rp, rs);
[b,a]=butter(n,wn,'low');
[h,w]=freqz(b,a,1024,fs);
y.Butterworth=filter(b,a,y.y);
%%
figure
plot(t,y.y,t,y1,t,y.Fourier,t,y.Butterworth);
xlabel('t')
ylabel('y')
legend('y','y1','Fourier','Butterworth')
grid on;grid minor;
saveas(gcf, '2 compare', 'png')

FourierFil.m

function filZ = FourierFil(fmax,fs,z)
num=length(z);
dt=1/fs;

%消除端点间段性（0点连续）
g0=zeros(1,num);
for(i=1:num)
    g0(i)=z(i)-z(1)-i*(z(num)-z(1))/(num-1);
end

%变为奇函数，拓展到双边 变成2*num-1个数
g=[-fliplr(g0) g0(2:num)];
         
%计算正弦级数系数
b=zeros(1,num-1);%第一个频率为0，故不需要计算第一个系数
T=zeros(1,num-1);
f=zeros(1,num-1);
for(k=1:num-1)%g()
    sumB=0;
    for(i=1:num-2) %g(1)=g(num)=0
        sumB=sumB+g0(i+1)*sin(k*pi*(i/(num-1)));%g(i)从g(1)开始        
    end
    b(k)=2/(num-1)*sumB;
    T(k)=2*(num-1)*dt/k;
    f(k)=1/T(k);    
end
figure
plot(f,abs(b))
xlabel('频率[Hz]')
ylabel('|b(k)|')

%反傅里叶变换
g1=zeros(1,num);
for(i=0:num-1)
    for(k=1:num-1)
        if(f(k)>fmax)%将大于截止频率的部分去除
            b(k)=0;
        end
        g1(i+1)=g1(i+1)+b(k)*sin(k*pi*(i/(num-1)));
    end    
end

figure
plot([1:num],g0,[1:num],g1);
legend('消除了端点间段性','消除了端点间段性傅里叶反变换')

%复原时间历程
for(i=0:num-1)
    filZ(i+1)=g1(i+1)+z(1)+i*(z(num)-z(1))/(num-1);
end

end