function y = cs_threshold(x, N, fs, M, fc, BW)
%%%%% 參量說(shuō)明 %%%%%
% x -- 信號(hào)
% N -- 循環(huán)譜檢測(cè)采樣長(zhǎng)度，必須小于等于信號(hào)序列長(zhǎng)度
% fs -- 采樣頻率, 檢測(cè)帶寬為-fs/2至fs/2
% M -- 平滑點(diǎn)數(shù), 時(shí)間分辨率*頻率分辨率=M
% BW -- 信號(hào)帶寬
%%%%%%%%%%%%%%%%%%%

%%% 參數(shù) %%% 
win = 'hamming'; % 平滑窗類型

d_alpha = fs/N; % 1/時(shí)間分辨率=循環(huán)頻率分辨率
alpha = 0:d_alpha:fs; % 循環(huán)頻率, 分辨率=1/時(shí)間分辨率
a_len = length(alpha); % 循環(huán)頻率取樣個(gè)數(shù)

f_len = floor(N/M-1)+1; % 最大平滑窗個(gè)數(shù), 即頻率采樣個(gè)數(shù)
f = -(fs/2-d_alpha*floor(M/2)) + d_alpha*M*(0:f_len-1); % 頻率采樣點(diǎn)位置

S = zeros(a_len, f_len); % 初始相關(guān)功率譜
i = 1; 

fl = ceil((2*fc-BW)/d_alpha)+1; % 所需檢測(cè)a的下限序號(hào)
rt = floor((2*fc+BW)/d_alpha)+1; % 所需檢測(cè)a的上限序號(hào)
peak_a = fl:rt; % a的檢測(cè)范圍
peak_f = (1:ceil(BW/d_alpha/M))+floor(f_len/2-floor(BW/d_alpha/M)/2); % f的最大檢測(cè)范圍
noise_n = 0; % 噪聲點(diǎn)數(shù)初始

%%% 信號(hào)fft變換 %%%
X = fftshift(fft(x(1:N))); 
X = X';

%%% 遍歷循環(huán)頻率取值 %%%
for alfa = alpha   
    
    if i >= fl && i <= rt % 僅檢測(cè)所需a范圍
        
        interval_f_N = round(alfa/d_alpha); % 循環(huán)頻率所對(duì)應(yīng)的頻譜序列序號(hào)
        f_N = floor((N-interval_f_N-M)/M)+1; % 平滑窗的個(gè)數(shù)
    
        %%% 生成平滑窗函數(shù) %%%
        g = feval(win, M); 
        window_M = g(:, ones(f_N,1));
    
        %%% 頻域序列平滑模板 %%%
        t = 1:M*f_N;
        t = reshape(t, M, f_N);
    
        %%% 計(jì)算X1,X2 %%%
        X1 = X(t).*window_M;
        X2 = X(t+interval_f_N).*window_M; 

        %%% 計(jì)算譜相關(guān) %%%
        St = conj(X1).*X2;
        St = mean(St, 1); % 平滑平均
    
        noise_n = noise_n + length(St); % 累計(jì)噪聲點(diǎn)
        S(i, floor((f_len-f_N)/2)+(1:f_N)) = St/N; % 將結(jié)果平移至序列中央以便作圖
    end
    
    i = i+1; 
end
%%% 遍歷循環(huán)頻率取值結(jié)束 %%%
 
%%% 循環(huán)功率譜作圖 %%%
%mesh(f, alpha, abs(S)); 
%axis tight;
%title('BPSK-OFDM');
%xlabel('f'); ylabel('a');

peak = S(peak_a, peak_f); % 取出峰值區(qū)域
S(peak_a, peak_f) = 0; 
noise = S(peak_a,:); % 取出噪聲區(qū)域 

p = sort(abs(peak(:)),'descend'); % 排序峰值區(qū)域
peak_mean = mean(p(1:ceil(length(p)*0.1))); % 前10%點(diǎn)平均

noise_n = noise_n-length(peak(:)); % 實(shí)際噪聲點(diǎn)數(shù)
n = sort(abs(noise(:)),'descend'); % 排序噪聲區(qū)域
noise_mean = mean(n(1:noise_n)); % 噪聲平均

y = peak_mean/noise_mean; % 計(jì)算比值