数学形态学 源码

% 数学形态学去噪
double u;
u=linspace（00.11280）;
y0=sin（2*pi*50*u）;%50HZ的信号
w=100*pi;

R=rand（12）;
T1=0.1*min（R）;T2=0.1*max（R）;
g1=heaviside（u-T1）;
g2=heaviside（u-T2）;
y=（1-0.8*（g1-g2））.*y0;  %该信号可作为原始信号

figure（12）
subplot（311）;
plot（y）;box off;xlabel（‘采样点‘）;title（‘理想扰动信号‘）

N=1280;%与信号同长度
noise=randn（1N）;
noise_mean=mean（noise）;%求均值
power_noise=var（noise）;
a1=sqrt（0.01/power_noise）;
noise=noise*a1;%将信号功率变为0.01

f=y+noise; %应该叠加噪声

% figure（1）
subplot（312）;
plot（f‘k‘）;box off;xlabel（‘采样点‘）;title（‘含噪声信号‘）

% g=zeros（1100）;
% for kk=1:100
%     g（kk）=sqrt（0.01-（（100-kk）/1000）^2）;
% end
% g
g=[00.04360.06000.07140.08000.08660.09170.09540.09800.09950.1000];
% g=0.05*ones（111）;
% g1=fliplr（g）;
% g=0.5*（g+g1）;
% kkk=sum（g）/length（g）;
% g=kkk*ones（111）;
% g=0.8*g;
N=1280;
M=length（g）;

%f关于g的腐蚀 f-g
fmg=zeros（size（f））;
fg=zeros（size（g））;
for i=1:N
    for j=1:M
        if i+j<=N
            
            fg（j）=f（i+j）-g（j）;
        end
    end
    fmg（i）=min（fg）;
end

%f关于g的膨胀 f+g
fpg=zeros（size（f））;
fg1=zeros（size（g））;
for i=1:N
    for j=1:M
        if i-j>=0
            if i-j>0
            
            fg1（j）=f（i-j）+g（j）;
            else
                fg1（j）=f（1）+g（j）;
            end
        end
    end
    fpg（i）=max（fg1）;
end

% f关于g的开运算f-g+g=fmg+g
fpg1=zeros（size（f））;
fg11=zeros（size（g））;
for i=1:N
    for j=1:M
        if i-j>=0
            if i-j>0
            
            fg11（j）=fmg（i-j）+g（j）;
            else
                fg11（j）=fmg（1）+g（j）;
            end
                
        end
    end
    fpg1（i）=max（fg11）;
end

%f+g-g
fmg1=zeros（size（f））;
fg12=zeros（size（g））;
for i=1:N
    for j=1:M
        if i+j<=N
            
            fg12（j）=fpg（i+j）-g（j）;
        end
    end
    fmg1（i）=min（fg12）;
end

%开运算fog.g=（f-g+g）