function bat1%非正交的二次樣條二進小波h = [0 0 0.125 0.375 0.375 0.125 0 0].*sqrt(2);g = [0 0 0 -1 1 0 0 0].*sqrt(2);g1 = [0 0.015625 0.109375 0.34375 -0.34375 -0.109375 -0.015625 0].*sqrt(2);%正交的二次樣條二進小波% h = [0 0 0 0 0.375 0.125 0.373].*sqrt(2);% g = [0 0 0 0 0.5798 0.0869 0.0061].*sqrt(2);%g1 = [0 0.015625 0.109375 0.34375 -0.34375 -0.109375 -0.015625 0].*sqrt(2);%正交的三次樣條二進小波% h = [0 0.0625 0.2500 0.3750 0.2500 0.0625 0].*sqrt(2);% g = [-0.00008 -0.01643 -0.10872 -0.59261 0 0.59261 0.10872 0.01643 0.00008].*sqrt(2);%g1 = [0 0.015625 0.109375 0.34375 -0.34375 -0.109375 -0.015625 0].*sqrt(2);l = 512;x = [1:l];n = 7;sa = 0;sn = 5;dsa = 2.5;snum = 12;noise = rand(1,l)-0.5;so(1:l/4) = 2;so(l/4+1:l*1/2) = 5;so(l*1/2+1:l*3/4) = 1;so(l*3/4+1:l) = 4;for i=1:snum    % s0 = sin(2*pi*x/l*4).*sin(2*pi*x/l*8).*sin(2*pi*x/l)*sa+2*sa;% s = s0+(rand(1,l)-0.5)*sn;% signalName = 's1';sa = sa + dsa;% s = sin(2*pi*x/l*2);s0 = so*sa;s = s0 + noise*sn;SignalEnergy(i) = sum(s0.*s0);NoiseEnergy(i) = sum((s-s0).*(s-s0));SNR1(i) = 10*log10(SignalEnergy(i)/NoiseEnergy(i));signalName = ['s3_SNR' num2str(SNR1(i))];noise1 = s-s0;noise1(1:l-1) = noise1(2:l)-noise1(1:l-1);      %求倒數noise1(1:l-1) =  noise1(2:l).*noise1(1:l-1);     %計算過零點NoiseVar(i) = length(find(noise1<0))/l;% s(1:l/4) = [1:l/4]./20 + l/8/20;% s(l/4+1:l*1/2) = s(l/4:-1:1);% s(l/2+1:l) = s(1:l/2);% s0 = s;% s = s0 + (rand(1,l)-0.5)*6;% signalName = 's2';picDir = ['pic' filesep signalName];if(exist(picDir,'dir'))    ii = 1;    picDir1 = [picDir '_' num2str(ii)];    while(exist(picDir1,'dir'))        ii = ii+1;        picDir1 = [picDir '_' num2str(ii)];    end    picDir = picDir1;endmkdir(picDir);% sr = ss(:,n);sr = SSNF(s,n,h,g,g1,picDir);NoiseEnergyAfter(i) = sum((sr'-s0).*(sr'-s0));SNR2(i) = 10*log10(SignalEnergy(i)/NoiseEnergyAfter(i));errorReduce(i) = (NoiseEnergy(i)-NoiseEnergyAfter(i))/NoiseEnergy(i);SNRadvanced(i) = SNR2(i) - SNR1(i);SNRrate(i) = (SNR2(i) - SNR1(i))/SNR1(i);noise2= sr'-s0;noise2(1:l-1) = noise2(2:l)-noise2(1:l-1);      %求倒數noise2(1:l-1) =  noise2(2:l).*noise2(1:l-1);     %計算過零點NoiseVarAfter(i) = length(find(noise2<0))/l;varReduce(i) = NoiseVar(i) - NoiseVarAfter(i);figuresubplot(6,1,1), hold on, title(['原始信號 E:' num2str(SignalEnergy(i))]),plot(s0), set(gca,'XTick',0:50:l);subplot(6,1,2), hold on, title(['噪聲信號 E:' num2str(NoiseEnergy(i))]),plot(s-s0), set(gca,'XTick',0:50:l);subplot(6,1,3), hold on, title(['含噪信號 SNR:'  num2str(SNR1(i))]),plot(s), set(gca,'XTick',0:50:l);subplot(6,1,4), hold on, title(['去噪信號 SNR:'  num2str(SNR2(i)) ' SNRadv: ' ...    num2str(SNRadvanced(i)) ' SNRrate: '  num2str(SNRrate(i))]),plot(sr), set(gca,'XTick',0:50:l);subplot(6,1,5), hold on, title(['誤差信號 E:'  num2str(NoiseEnergyAfter(i)) ' error redue: ' ...    num2str(errorReduce(i))]),plot(sr-s0'), set(gca,'XTick',0:50:l);subplot(6,1,6), hold on, title(['誤差變化 E:'  num2str(NoiseVarAfter(i)) ' var redue: ' ...    num2str(varReduce(i))]),plot(noise2), set(gca,'XTick',0:50:l);saveas(gcf, [picDir filesep 'finalResult'], 'png'), close(gcf);endx = [1:snum];figurehold on;[AX,H1,H2] = plotyy(x,SNR1,x,errorReduce,'plot');set(get(AX(1),'Ylabel'),'String','SNR')set(get(AX(2),'Ylabel'),'String','ErrorReduced')set(H1,'LineStyle','--')set(H2,'LineStyle','-.')title('去噪效果')saveas(gcf, ['pic' filesep 'SNRvsError'], 'png')%, close(gcf);figurehold on;[AX,H1,H2] = plotyy(x,SNR1,x,varReduce,'plot');set(get(AX(1),'Ylabel'),'String','SNR')set(get(AX(2),'Ylabel'),'String','VarReduced')set(H1,'LineStyle','--')set(H2,'LineStyle','-.')title('去噪效果')saveas(gcf, ['pic' filesep 'SNRvsVarReduce'], 'png')%, close(gcf);figurehold on;plot(x,NoiseVar,'-',x,NoiseVarAfter,'-.',x,varReduce,'--');legend('VAR orginal','VAR after denoising','VAR reduced','Location','East');title('VAR比較')saveas(gcf, ['pic' filesep 'Varcompare'], 'png')figurehold on;plot(x,SNR1,'-',x,SNR2,'-.',x,SNRadvanced,'--');legend('SNR orginal','SNR after denoising','SNR advanced','Location','East');title('SNR比較')saveas(gcf, ['pic' filesep 'SNRcompare'], 'png')