参考文献
[1]刘冲,马立修,潘金凤,等.联合VMD与改进小波阈值的局放信号去噪[J].现代电子技术,2021,44(21):45
参考文献引言
改进原理(三)
部分核心代码
% 测试文件clc;clear;close all%% 1.生成仿真信号N=1000;Fs=1000;%%采样频率自己设置t=((0:N-1)*1/Fs)';f1=5;f2=30;y1=5*sin(2*pi*f1*t);y2=3*sin(2*pi*f2*t);y=y1+y2;SNR=10;rng(100);sig_noise = awgn(y,SNR,'measured'); %添加噪声%% 画图ylabel('\fontname{宋体}幅值');xlabel('\fontname{Times new roman}Frequency/\it{Hz}');subplot(221);plot(y,'k');title('\fontname{宋体}原始信号');ylabel('\fontname{宋体}幅值');xlabel('采样点数/n');subplot(222);pFFT(y,Fs);title('\fontname{宋体}频谱图');ylabel('\fontname{宋体}幅值');xlabel('\fontname{Times new roman}Frequency/\it{Hz}');subplot(223);plot(sig_noise,'k');title('\fontname{宋体}加噪后信号');ylabel('\fontname{宋体}幅值');xlabel('采样点数/n');subplot(224);pFFT(sig_noise,Fs);title('\fontname{宋体}频谱图');ylabel('\fontname{宋体}幅值');xlabel('\fontname{Times new roman}Frequency/\it{Hz}');%% 2.小波阈值去噪%%小波基函数的选择 %%分解层数的选择 %%阈值选择规则 %%小波基函数的选择：https://www.mathworks.cn/help/wavelet/ref/wfilters.html?searchHighlight=wname&s_tid=srchtitle_wname_2#d123e130597wname='db3'; lev=4; tptr='heursure'; %'rigsure','heursure','sqtwolog','minimaxi'option=[];SORH={'h','s','a2'};for i=1:length(SORH)WT_sig{i}=IWTD(sig_noise,wname,SORH{i},lev,tptr,option); %调用函数进行滤波endfigure('color','w');set(gcf, 'Position', [400 300 500 500]);subplot(2,2,1);plot(sig_noise,'k'); title('\fontname{宋体}加噪后信号');ylabel('幅值');xlabel('采样点数/n');Name={'硬阈值小波去噪','软阈值小波去噪','改进阈值函数小波去噪'};for i=1:length(SORH)subplot(2,2,i+1);plot(WT_sig{i});ylabel('\fontname{宋体}幅值');xlabel('采样点数/n');title(Name{i});endSimulate_Signal=y;for i=1:length(SORH)[~,RMSE,RMS,MAE,SNR,~,cross_core,R,SER,NM,~,BIAS]=evaluate_denoising_metrics(Simulate_Signal,WT_sig{i});disp(['---------------------------------',SORH{i},'小波的去噪的指标----------------------------------']);disp(['均方根误差RMSE=',num2str(RMSE),'均方根值RMS=',num2str(RMS),'，平均绝对值误差MAE=',num2str(MAE),'，信噪比SNR=',num2str(SNR),'，互相关系数cc=',num2str(cross_core),'平滑度指标R=',num2str(R),'，信号能量比SER=',num2str(SER),'，噪声的模NM=',num2str(NM),'，信号偏差BIAS=',num2str(BIAS)]);endfunction thr = kthselect(x,tptr,lev,wname)validateattributes(x, {'numeric'}, {'2d'}, 'THSELECT', 'X')if isStringScalar(tptr)tptr = convertStringsToChars(tptr);endif isrow(x)x = x(:);end[n,m] = size(x);switch tptrcase 'rigrsure'sx = sort(abs(x),1);sx2 = sx.^2;N1 = repmat((n-2*(1:n))',1,m);N2 = repmat((n-1:-1:0)',1,m);CS1 = cumsum(sx2,1);risks = (N1+CS1+N2.*sx2)./n;[~,best] = min(risks,[],1, 'linear');% thr will be row vectorthr = sx(best);case 'heursure'%%hthr = (2*log(n)).^0.5;eta = sum(abs(x).^2-n,1)./n;crit = (log(n)/log(2))^(1.5)/(n.^0.5);thr = thselect(x,'rigrsure');thr(thr > hthr) = hthr;thr(eta < crit) = hthr;case 'sqtwolog'thr = (2*log(n)).^0.5;thr = repelem(thr,size(x,2));case 'minimaxi'if n <= 32thr = repelem(0, m);elset = 0.3936 + 0.1829*(log(n)/log(2));thr = repelem(t,m);endcase 'visushrink'[c,l]=wavedec(x,lev,wname);for i = 1:levd{i} = detcoef(c,l,i); % 提取第i层的高频系数thr(i) = median(abs(d{i}))/0.6745*sqrt(2*log(length(d{i})));endotherwiseerror(message('Wavelet:FunctionArgVal:Invalid_ArgVal'))endend
