%%初始化程序
clear,clc
t1=clock;
%% 载入7.6日00:15到7.21日24:00的数据
load('SJ76721.mat');%matrix
YSJ= SJ76721;
disp('------7.6日00:15到7.21日24:00的降噪处理结果-------');
% 载入7.6日00:15到7.24日24:00的数据
% load('SJ76724.mat');
% YSJ= SJ76724;
% disp('------7.6日00:15到7.24日24:00的降噪处理结果-------');
%% 数据预处理
[c,l]=size(YSJ);
Y=[];
for i=1:c
Y=[Y,YSJ(i,:)];
end
[c1,l1]=size(Y);
X=[1:l1];
%% 绘制原始信号图像
figure(1);
plot(X,Y);
xlabel('负荷个数/(个)');
ylabel('kW·h');
title('原始信号');
%% 硬阈值处理
lev=3;
xd=wden(Y,'heursure','h','one',lev,'db4');%硬阈值去噪处理后的信号序列
figure(2)
plot(X,xd)
xlabel('负荷个数/(个)');
ylabel('kW·h');
title('硬阈值去噪处理')
set(gcf,'Color',[1 1 1])
%% 软阈值处理
lev=3;
xs=wden(Y,'heursure','s','one',lev,'db4');%软阈值去噪处理后的信号序列
figure(3)
plot(X,xs)
xlabel('负荷个数/(个)');
ylabel('kW·h');
title('软阈值去噪处理')
set(gcf,'Color',[1 1 1])
%% 改进阈值后的去噪处理
lev=3;
xz=wden(Y,'sqtwolog','s','sln',lev,'db4');%改进阈值去噪处理后的信号序列
figure(4)
plot(X,xz);
xlabel('负荷个数/(个)');
ylabel('kW·h');
title('改进阈值后的去噪处理')
set(gcf,'Color',[1 1 1])
%% 计算信噪比SNR
Psig=sum(Y*Y')/l1;
Pnoi1=sum((Y-xd)*(Y-xd)')/l1;
Pnoi2=sum((Y-xs)*(Y-xs)')/l1;
Pnoi3=sum((Y-xz)*(Y-xz)')/l1;
SNR1=10*log10(Psig/Pnoi1);
SNR2=10*log10(Psig/Pnoi2);
SNR3=10*log10(Psig/Pnoi3);
%% 计算均方根误差RMSE
RMSE1=sqrt(Pnoi1);
RMSE2=sqrt(Pnoi2);
RMSE3=sqrt(Pnoi3);
%% 计算均方根误差RMSE
disp(['硬阈值去噪处理的SNR=',num2str(SNR1),',RMSE=',num2str(RMSE1)]);
disp(['软阈值去噪处理的SNR=',num2str(SNR2),',RMSE=',num2str(RMSE2)]);
disp(['改进阈值后的去噪处理SNR=',num2str(SNR3),',RMSE=',num2str(RMSE3)]);
t2=clock;
tim=etime(t2,t1);
disp(['------------------运行耗时',num2str(tim),'秒-------------------'])