19、matlab信号预处理中的中值滤波(medfilt1()函数)和萨维茨基-戈雷滤波滤(sgolayfilt()函数)

1、中值滤波：medfilt1()函数

说明：一维中值滤波

1）语法

语法1：y = medfilt1(x) 将输入向量x应用3阶一维中值滤波器。

语法2：y = medfilt1(x,n) 将一个n阶一维中值滤波器应用于x。

语法3：y = medfilt1(x,n,[],dim) 指定过滤器操作的维度dim。

2）参数说明

x:输入信号 y:输出信号 dim:筛选维度

3）NaN(信号数据缺失)处理

'includenan' ：返回过滤信号，以便包含NaN的任何段的中值也是NaN。

'omitnan' 返回过滤后的信号，使得包含 NaN 的任何段的中值为非 NaN 值的中值。如果一个段的
所有元素都是 NaN，则结果为 NaN终点过滤，指定为'零填充'或'截断'。

'zeropad' 在端点之外，信号被视为零。 'truncate' 在接近信号边缘时计算较小段的中位数。

2、中值滤波实验

1）通过中值滤波进行降噪

代码

fs = 200;%频率
t = 0:1/fs:2;
x = sin(2*pi*t*5)+0.25*sin(2*pi*t*80);
y = medfilt1(x,3);
plot(t,x,'-','color','r')
hold on;
plot(t,y,'-^','color','g')
legend('原始信号','滤波后信号')

视图效果

2）带有尖峰和丢失样本的多通道信号中值滤波

说明：

生成一个由不同频率正弦波组成的双通道信号。在随机位置加入尖峰。在随机位置用NaN添加缺失样本。重置随机数生成器设置加入噪声位置。

多通道信号生成代码

rng('default')
n = 59;
x = sin(pi./[10 20]'*(1:n)+pi/6)';
spk = randi(100,9,1);
x(spk) = x(spk)*2;
x(randi(100,6,1)) = NaN;
plot(x,'color','r')
legend('信号1','信号2')

试图效果

中值滤波效果对比

代码

rng('default')
n = 59;
x = sin(pi./[10 20]'*(1:n)+pi/6)';
spk = randi(100,9,1);
x(spk) = x(spk)*2;
x(randi(100,6,1)) = NaN;
plot(x,'color','r')
legend('信号1','信号2')
hold on;
y = medfilt1(x,8);
plot(y,'color','g')
legend('信号1中值滤波','信号2中值滤波')

试图效果

信号缺失部分处理（NaN）

代码

y = medfilt1(x,4,'omitnan');
plot(y)

边缘滤波

代码

y = medfilt1(x,4,'omitnan','truncate');
plot(y)

视图效果

3、萨维茨基-戈雷滤波滤波器：sgolayfilt()

说明：萨维茨基-戈雷滤波

语法

语法1：y = sgolayfilt(x,order,framelen) 对向量 x 中的数据应用多项式阶数为 order、帧长度为
framelen 的萨维茨基-戈雷有限冲激响应 (FIR) 平滑滤波器。
语法2：y = sgolayfilt(x,order,framelen,weights) 指定在最小二乘最小化过程中要使用的加权向量。
语法3：y = sgolayfilt(x,order,framelen,weights,dim) 指定滤波器沿其运算的维度。

参数

x：输入信号 order：多项式阶数 framelen：帧长度 weights：加权数组 dim:要沿其滤波的维度

1）稳态和瞬变萨维茨基-戈雷滤波器

萨维茨基-戈雷滤波器滤波代码

order = 4;%参数设置 可以根据需求设置
framelen = 13;
l = 40;
x = randn(l,1);
sgf = sgolayfilt(x,order,framelen);
plot(x,':')
hold on
plot(sgf,'.-')
legend('原信号','戈雷滤波信号')
hold on

视图效果

2）稳态萨维茨基-戈雷滤波器滤波代码

order = 4;%参数设置 可以根据需求设置
framelen = 13;
l = 40;
x = randn(l,1);
sgf = sgolayfilt(x,order,framelen);
plot(x,':')
hold on
plot(sgf,'.-')
% legend('原信号','戈雷滤波信号')
hold on
m = (framelen-1)/2;
B = sgolay(order,framelen);
steady = conv(x,B(m+1,:),'same');%原信号与过滤信号卷积滤波 得到稳态部分
plot(steady)
legend('原信号','戈雷滤波信号','稳态部分')

视图效果

3）启动瞬态和终止瞬态代码

order = 4;%参数设置 可以根据需求设置
framelen = 13;
l = 40;
x = randn(l,1);
sgf = sgolayfilt(x,order,framelen);
plot(x,':')
hold on
plot(sgf,'.-')
% legend('原信号','戈雷滤波信号')
hold on
m = (framelen-1)/2;
B = sgolay(order,framelen);
steady = conv(x,B(m+1,:),'same');%原信号与过滤信号卷积滤波 得到稳态部分
plot(steady)
legend('原信号','戈雷滤波信号','稳态部分')

ybeg = B(1:m,:)*x(1:framelen);%启动瞬态
yend = B(framelen-m+1:framelen,:)*x(l-framelen+1:l);%终止瞬态
cmplt = steady;
cmplt(1:m) = ybeg;
cmplt(l-m+1:l) = yend;

plot(cmplt)
legend('原信号','戈雷滤波信号','稳态部分','完整信号')