This article addresses the problem of estimating the spectral correlation function (SCF), which provides quantitative characterization in the frequency domain of wide-sense cyclostationary properties of random processes which are considered to be the theoretical models of observed time series or discrete-time signals. The theoretical framework behind the SCF estimation is briefly reviewed so that an important difference between the width of the resolution cell in bifrequency plane and the step between the centers of neighboring cells is highlighted. The outline of the proposed double-number fast Fourier transform algorithm (2N-FFT) is described in the paper as a sequence of steps directly leading to a digital signal processing technique. The 2N-FFT algorithm is derived from the time-smoothing approach to cyclic periodogram estimation where the spectral interpolation based on doubling the FFT base is employed. This guarantees that no cyclic frequency is left out of the coverage grid so that at least one resolution element intersects it. A numerical simulation involving two processes, a harmonic amplitude modulated by stationary noise and a binary-pulse amplitude-modulated train, demonstrated that their cyclic frequencies are estimated with a high accuracy, reaching the size of step between resolution cells. In addition, the SCF components estimated by the proposed algorithm are shown to be similar to the curves provided by the theoretical models of the observed processes. The comparison between the proposed algorithm and the well-known FFT accumulation method in terms of computational complexity and required memory size reveals the cases where the 2N-FFT algorithm offers a reasonable trade-off.