传统数学形态学滤波器与广义数学形态学滤波器的matlab程序实现.rar

%% 2019年3月14日，测试广义数学形态学效果主程序-脉冲 %这里的腐蚀、膨胀结构元素长度不一致 %对于脉冲噪声，单一的传统开、闭，效果不好，开可以滤除正脉冲，闭可以滤除负脉冲 %但是对于广义的开、闭，效果改善了，广义开的上升部分，比原始信号右移了L1-1个单位；广义闭的下降部分，比原始信号左移了L1-1个单位 %对于OC\CO\GCO\GOC效果都很好，接近原始波形，区别不是很大。 %% 0-清空工作空间 clear close all clc %% 1-建立原始波形 load('sinx.mat'); %% 2-加噪声 %正脉冲噪声 x=sinx(1:5000);%全局的输入函数，原始信号 am=0.5;%脉冲幅值为0.5v wm=2;%脉冲宽度设置为2个采样点 %加入正脉冲噪声 y=pulsenoise(x,501,am,wm); y=pulsenoise(x,751,am,wm); y=pulsenoise(y,1001,am,wm); y=pulsenoise(y,1251,am,wm); y=pulsenoise(y,1501,am,wm); y=pulsenoise(y,1751,am,wm); y=pulsenoise(y,2001,am,wm); %加入负脉冲噪声 y=pulsenoise(y,2501,-am,wm); y=pulsenoise(y,2751,-am,wm); y=pulsenoise(y,3001,-am,wm); y=pulsenoise(y,3251,-am,wm); y=pulsenoise(y,3501,-am,wm); y=pulsenoise(y,3751,-am,wm); y=pulsenoise(y,4001,-am,wm); %查看效果 figure t=[0.000001:0.000001:0.000001*length(x)]; plot(t,y,'k'); xlabel('时间t/s'); ylabel('幅值/V'); title('加噪声后波形') %% 3-设置参数 %本程序的目的是找到最佳的结构元素形状和结构宽度，里面参数均可以修改 ORIA=0.0; %A的起点******* ENDA=0.0; %A的终点****** STEPA=0.01; %A的步长********* A=[ORIA:STEPA:ENDA]; %A矩阵，幅值矩阵 aa=length(A); %aa表示矩阵的列，和A的宽度一样大 ORIL1=2; %L1的起点，2开始************* ENDL1=20; %L1的终点********** ORIL2=2*ORIL1-1; ENDL2=2*ENDL1-1; %L2的终点********** L1=[ORIL1:1:ENDL1]; %L1步长固定为1 L2=[ORIL2:2:ENDL2]; %L2步长固定为2 bb1=length(L1); %bb表示矩阵的行，和L的宽度一样大 bb2=length(L2); %bb表示矩阵的行，和L的宽度一样大 OMAE=zeros(aa,bb1); CMAE=zeros(aa,bb1); GCMAE=zeros(aa,bb1); GOMAE=zeros(aa,bb1); OCMAE=zeros(aa,bb1); COMAE=zeros(aa,bb1); GOCMAE=zeros(aa,bb1); GCOMAE=zeros(aa,bb1); YMAE=zeros(aa,bb1); GYMAE=zeros(aa,bb1); OSNR=zeros(aa,bb1); CSNR=zeros(aa,bb1); GCSNR=zeros(aa,bb1); GOSNR=zeros(aa,bb1); OCSNR=zeros(aa,bb1); COSNR=zeros(aa,bb1); GOCSNR=zeros(aa,bb1); GCOSNR=zeros(aa,bb1); YSNR=zeros(aa,bb1); GYSNR=zeros(aa,bb1); %% 结构元素 for k=1:1:aa for j=1:1:bb1 %直线 for i=1:L1(j) b1(i)=A(k);%将当前A的值赋给b1 end for i=1:L2(j) b2(i)=A(k);%将当前A的值赋给b2 end % plot(b1); % plot(b2); %余弦 % cos_x=(-pi/2):pi/(L1(j)-1):(pi/2); % cos_y=A(k)*cos(cos_x); % b1=roundn(cos_y,-4);%取小数点后四位，便于运算 % % cos_x1=(-pi/2):pi/(L2(j)-1):(pi/2); % cos_y1=A(k)*cos(cos_x1); % b2=roundn(cos_y1,-4);%取小数点后四位，便于运算 % plot(b1); % plot(b2); %半圆 % cir_x=-A(k):(2.*A(k)./(L1(j)-1)):A(k); % cir_y=sqrt(A(k)^2-cir_x.^2); % b1=roundn(cir_y,-4);%取小数点后四位，便于运算 % % cir_x1=-A(k):(2.*A(k)./(L2(j)-1)):A(k); % cir_y1=sqrt(A(k)^2-cir_x1.^2); % b2=roundn(cir_y1,-4);%取小数点后四位，便于运算 % plot(b1); % plot(b2); %% 三角1(2N时候适用) % 三角(L为奇数) % if mod(L1(j),2)==1 % mid=(L1(j)-1)/2; % tt=0:1:mid; % b1=A(k)./mid.*tt; % for i=1:mid % b1(mid+1+i)=b1(mid+1-i); % end % % mid1=(L2(j)-1)/2; % tt1=0:1:mid1; % b2=A(k)./mid1.*tt1; % for i=1:mid1 % b2(mid1+1+i)=b2(mid1+1-i); % end % else % %三角(L为偶数) % mid=(L1(j)-1)/2; % po=round(mid); % tt=0:1:fix(mid);%向零取整 % b1=A(k)/mid.*tt; % for i=1:po % b1(po+i)=b1(po+1-i); % end % % mid1=(L2(j)-1)/2; % po1=round(mid1); % tt1=0:1:fix(mid1);%向零取整 % b2=A(k)/mid1.*tt1; % for i=1:po1 % b2(po1+i)=b2(po1+1-i); % end % end % plot(b1); % plot(b2); %% 三角1(2N-1时候适用) % 三角(b1为奇数) % if mod(L1(j),2)==1 % mid=(L1(j)-1)/2; % tt=0:1:mid; % b1=A(k)./mid.*tt; % for i=1:mid % b1(mid+1+i)=b1(mid+1-i); % end % %b2为偶数 % mid1=(L2(j)-1)/2; % po1=round(mid1); % tt1=0:1:fix(mid1);%向零取整 % b2=A(k)/mid1.*tt1; % for i=1:po1 % b2(po1+i)=b2(po1+1-i); % end % else % %三角(b1为偶数) % mid=(L1(j)-1)/2; % po=round(mid); % tt=0:1:fix(mid);%向零取整 % b1=A(k)/mid.*tt; % for i=1:po % b1(po+i)=b1(po+1-i); % end % %b2为奇数 % mid1=(L2(j)-1)/2; % tt1=0:1:mid1; % b2=A(k)./mid1.*tt1; % for i=1:mid1 % b2(mid1+1+i)=b2(mid1+1-i); % end % end % plot(b1); % plot(b2); %% 4-测试 insignal=y; O2=fushi2(insignal,b1,0); O2=pengzhang2(O2,b1,0); %移动波形（掐头去尾） sinx0=x; sinx0(1:(L1(j)-1))=[];%掐头 O3=O2; len=length(O3); O3((len+1-(L1(j)-1)):len)=[];%去尾 % figure % plot(sinx0,'k'); % hold on % plot(O3,'r'); % hold off % legend('原始信号','O运算'); % title('O效果') C2=pengzhang2(insignal,b1,0); C2=fushi2(C2,b1,0); sinx0=x; sinx0(1:(L1(j)-1))=[];%掐头 C3=C2; len=length(C3); C3((len+1-(L1(j)-1)):len)=[];%去尾 % figure % t1=[0.000001:0.000001:0.000001*length(sinx0)]; % plot(t1,sinx0,'k'); % hold on % plot(t1,C3,'r'); % hold off % xlabel('时间t/s'); % ylabel('幅值/V'); % legend('原始信号','C运算'); % title('C效果') GO2=GO(insignal,b1,b2); GC2=GC(insignal,b1,b2); sinx0=x; sinx0(1:(L1(j)-1))=[];%掐头 GO3=GO2; len=length(GO3); GO3((len+1-(L1(j)-1)):len)=[];%去尾 GC3=GC2; len=length(GC3); GC3((len+1-(L1(j)-1)):len)=[];%去尾 % figure % plot(sinx0,'k'); % hold on % plot(GO3,'r'); % hold off % legend('原始信号','GO运算'); % title('GO效果') % % figure % t1=[0.000001:0.000001:0.000001*length(sinx0)]; % plot(t1,sinx0,'k'); % hold on % plot(t1,GC3,'r'); % hold off % xlabel('时间t/s'); % ylabel('幅值/V'); % legend('原始信号','GC运算'); % title('GC效果') insignal=y; OC2=OC(insignal,b1); %移动波形（掐头去尾） sinx1=x; sinx1(1:2*(L1(j)-1))=[];%掐头 OC3=OC2; len=length(OC3); OC3((len+1-2*(L1(j)-1)):len)=[];%去尾 % figure % plot(sinx1,'k'); % hold on % plot(OC3,'r'); % hold off % legend('原始信号','OC运算'); % title('OC效果') % insignal=y; GOC2=GOC(insignal,b1,b2); %移动波形（掐头去尾） sinx2=x; sinx2(1:3*(L1(j)-1))=[];%掐头 GOC3=GOC2; len=length(GOC3); GOC3((len+1-3*(L1(j)-1)):len)=[];%去尾 % figure % plot(sinx2,'k'); % hold on % plot(GOC3,'r'); % hold off % legend('原始信号','GOC运算'); % title('GOC效果') insignal=y; CO2=CO(insignal,b1); %移动波形（掐头去尾） sinx1=x; sinx1(1:2*(L1(j)-1))=[];%掐头 CO3=CO2; len=length(CO3); CO3((len+1-2*(L1(j)-1)):len)=[];%去尾 % figure % plot(sinx1,'k'); % hold on % plot(C