当前位置：首页 > news >正文

【小波去噪】【matlab】基于小波分析的一维信号滤波（对照组：中值滤波、均值滤波、高斯滤波）

news 来源：原创 2024/9/20 11:44:22

链接1-傅里叶变换
链接2-傅立叶分析和小波分析间的关系
链接3-小波变换（wavelet transform）的通俗解释
链接4-小波基的选择

1.示例代码

function main_wavelet
clc
clear
close all
warning off
%% 1.信号生成
time_length = 10;%总时长，秒
sample_rate = 100;% 采样率，Hz
freq = 10; % 生成的信号的频率，Hz
t = linspace(0, time_length, time_length * sample_rate);  % 生成时间轴
signal = 2*sin(2 * pi * freq * t);% 生成正弦信号
noise = 0.5 * randn(size(signal));% 生成白噪声，高斯过程
noisy_signal = signal + noise;% 将正弦信号和噪声信号相加生成含噪采集数据figure(1)
plot(t,signal,'b-')
hold on
plot(t, noisy_signal,'r-');
hold on
plot(t, noise,'k-');
xlabel('Time (s)');
ylabel('Amplitude');
legend('原信号','含噪信号','噪声')% save('noise.mat','noisy_signal')
%% 2.小波变换
% 设置小波函数和变换阶数
wname = 'db4';  % 选用 Daubechies 4 小波
level = 2;      % 小波变换的阶数
% 进行小波变换
[C, L] = wavedec(noisy_signal, level, wname);
% 提取细节系数
D = detcoef(C, L, level); 
% 对细节系数进行阈值处理
sigma = median(abs(D)) / 0.6745;  % 计算阈值，robust estimator估计方法
D = wthresh(D, 'h', sigma);       % 硬阈值处理
% 重构信号
data_denoised = wrcoef('a', C, L, wname, level);%% 3.绘制结果
figure(3)
subplot(3,1,1); plot(t, noisy_signal,'b-'); title('原始信号(含噪声)');
subplot(3,1,2); plot(t,signal,'b-'); title('原始信号(无噪声)');
subplot(3,1,3); plot(t,data_denoised,'b-'); title('降噪后的信号');figure(4)
plot(t,signal,'b-')
hold on
plot(t,data_denoised,'r-')
legend('原始信号(无噪声)','降噪后的信号')

2.对比滤波器

作为对比，以下代码直接在第1节的main_wavelet.m之后运行即可
中值滤波、均值滤波、高斯滤波

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 一维中值滤波
n=21;                % n为模板长度，值可以改变
y=medfilt1(noisy_signal,n);
figure(101)
plot(t,noisy_signal,'b-')
hold on
plot(t,y,'r-')
legend('中值滤波前的序列','中值滤波后的序列')
%% 一维均值滤波
n=5;                 % n为模板长度，值可以改变
mean=ones(1,n)./n;   % mean为1×n的模板，各数组元素的值均为1/n
y = conv(noisy_signal,mean);
y=y(1:length(y)-length(mean)+1);
figure(102);
plot(t,noisy_signal,'b-')
hold on
plot(t,y,'r-')
legend('均值滤波前的序列','均值滤波后的序列')
%% 一维高斯滤波
gau=[0.0009 0.0175 0.1295 0.3521 0.3521 0.1295 0.0175 0.0009];   % 标准差为1时的高斯函数一维模板，如果标准差不为1，则要修改模板
% gau=[0.0090,0.4910,0.4910,0.0090];                               % 标准差为0.5时的高斯函数一维模板，如果标准差不为0.5，则要修改模板
y=conv(noisy_signal,gau);
y=y(1:length(y)-length(gau)+1);
figure(103);
plot(t,noisy_signal,'b-')
hold on
plot(t,y,'r-')
legend('高斯滤波前的序列','高斯滤波后的序列')

3.封装代码

主代码

function wavelet_fun_main240919
clc
clear
close all
warning off
%% 1.信号生成
time_length = 10;%总时长，秒
sample_rate = 100;% 采样率，Hz
freq = 10; % 生成的信号的频率，Hz
t = linspace(0, time_length, time_length * sample_rate);  % 生成时间轴
noise = 0.5 * randn(size(t));% 生成白噪声，高斯过程
Results_JGLD_ideal = 2*sin(2 * pi * freq * t);% 生成正弦信号
Results_JGLD_error = Results_JGLD_ideal + noise;% 将正弦信号和噪声信号相加生成含噪采集数据%% 小波降噪处理
dt=1/sample_rate;%采样时间间隔
data_denoised=WaveletGenerate(Results_JGLD_ideal,Results_JGLD_error,dt);%% 画图
figure(1)
plot(t,Results_JGLD_ideal,'b-')
hold on
plot(t,data_denoised,'r-')
legend('原始信号(无噪声)','降噪后的信号')

WaveletGenerate.m

function data_denoised=WaveletGenerate(Results_JGLD_ideal,Results_JGLD_error,dt)
%% 1.信号生成
signal = Results_JGLD_ideal;% 原始信号
noisy_signal = Results_JGLD_error;%含噪采集数据
% dt=0.01;%时间间隔
t=dt:dt:(length(signal))*dt;
%% 2.小波变换
% 设置小波函数和变换阶数
wname = 'db4';  % 选用 Daubechies 4 小波
level = 1;      % 小波变换的阶数
% 进行小波变换
[C, L] = wavedec(noisy_signal, level, wname);
% 提取细节系数
D = detcoef(C, L, level); 
% 对细节系数进行阈值处理
sigma = median(abs(D)) / 0.6745;  % 计算阈值，robust estimator估计方法
D = wthresh(D, 'h', sigma);       % 硬阈值处理
% 重构信号
data_denoised = wrcoef('a', C, L, wname, level);