基于大津算法的峰值-非峰值自适应阈值算法

基于大津算法的峰值-非峰值自适应阈值算法

function [ threshold ] = otus1D( signal )
%输入参数：一维信号序列
%输出： 阈值
    res = tabulate(signal(:));
    histogram = res(:,2)';
    keys = res(:,1)';
    threshold = 0;
    sum0 = 0;
    sum1 = 0;
    cnt0 = 0;
    cnt1 = 0;
    w0 = 0;
    w1 =0;
    u0 = 0;
    u1 = 0;
    variance = 0;
    u =0;
    maxVariance = 0;
    
    for i = keys(2:end)
            sum0 = 0;
            sum1 = 0;
            cnt0 = 0;
            cnt1 = 0;
            w0 = 0;
            w1 =0;
            
            for j =keys
                if(j<i)
                    index = find(j ==keys);
                    cnt0  = cnt0 + histogram(index);
                    sum0 = sum0 +j*histogram(index);
                end
            end
            u0 = sum0/cnt0;
            w0 = cnt0/length(signal);
             for k =keys
                if(k>=i)
                    index1 = find(k ==keys);
                    cnt1  = cnt1 + histogram(index1);
                    sum1= sum1+k*histogram(index1);
                end
             end
            u1 = sum1/cnt1;
            w1 = 1-w0;
            u = u0*w0 + u1*w1;
            variance = w0*(u0-u)^2 + w1*(u1-u)^2;
            if(variance > maxVariance)
                maxVariance = variance;
                threshold = i;
            end
    end
    

end

测试如下：

clear all;
Fs = 1000;
M = 1000;
T = 1/Fs;
t = (1:M-1).*T;

y = sin(2*pi*300.*t)+cos(2*pi*30.*t);
k = find(y>0);
positive_y = y(k);
plot(positive_y);
hold on
threshold = otus1D(positive_y);
plot([0,length(positive_y)-1],[threshold,threshold])

测试结果