R&B test

function rbtest

DEBUG = 0;

% number of samples per second

fs = [1000 100000];

% snr

snr = [1:3:120];

f_std = zeros(length(fs), length(snr));

% number of iterations

iter = 2;

% number of steps

nsteps = 25;

% initial zero padding

K = 2;

for w = [1:length(fs)]

    f_bias = zeros(1, iter);

    for q = [1:length(snr)]

        for i = [1:iter]

            % generate random frequency and phase

            fr = rand() * fs(w) / 2;
    %         fr = fs(w)/4;

            ph = rand() * 2*pi;

            % generate test signal

            ideal_sgn = sqrt(2)*sin([0:fs(w)-1]/fs(w) * 2*pi * fr + ph);

            sgn = awgn (ideal_sgn, snr(q) - 10*log10(fs(w)));

            fourier_abs = abs(fft([sgn zeros(1, length(sgn)*(K-1))]));
    %         fourier_abs = abs(fft([hann(length(sgn))'.*sgn zeros(1, length(sgn)*(K-1))]));
            [my, mi] = max(fourier_abs(1:end/2));
            xc = [mi-2 mi]/K;
            for st = [1:nsteps]
                yout = dichotomy(xc);
                if (yout(3) > yout(1))
                    xc = [(xc(2)+xc(1))/2 xc(2)];
                else
                    xc = [xc(1) (xc(2)+xc(1))/2];
                end
            end

            xres = (xc(2)+xc(1))/2;

            % taking possible overflow into account
            if (xres < 0)
                xres = -xres;
            elseif (xres >= fs(w)/2)
                xres = fs(w) - xres;
            end

            f_bias(i) = xres - fr;

            if (DEBUG == 1)
                fr
                mi
                plot([0:fs(w)*K-1]/K, fourier_abs);
            end

        end

        % sum(f_bias)/length(f_bias)
        f_std(w, q) = sqrt(sum(f_bias.^2)/length(f_bias));
    end
end

figure
    semilogy(snr, f_std(1, :), 'k-x', snr, f_std(2, :), 'b-*', snr, sqrt(3/(2*pi^2)*10.^(-0.1*snr)), 'r-o');
    xlabel ('SNR (1 sec), dB');
    ylabel ('Standard deviation, Hz')
    legend (sprintf('fs = %d Hz', fs(1)), sprintf('fs = %d Hz', fs(2)), 'MCRB');
    title (sprintf('Standard deviation of frequency estimation for a real tone of 1 sec duration,\n%d iterations, FFT with K=%d zero padding, %d dichotomy steps', iter, K, nsteps))


function yo = dichotomy(xi)
    if (length(xi) ~= 2)
        display('error function dichotomy');
        return;
    end

    xin = [xi(1) (xi(1)+xi(2))/2 xi(2)];

    yo = abs( sgn * exp(xin'*j*[0:fs(w)-1]/fs(w) * 2*pi  )' );
end

end