function?y=cs（xNfsM）
%%%%%?參量說明?%%%%%參量怎么取？
%?x?--?信號
%?N?--?循環譜檢測采樣長度，必須小于等于信號序列長度
%?fs?--?采樣頻率?檢測帶寬為-fs/2至fs/2
%?M?--?平滑點數?時間分辨率*頻率分辨率=M
%%%%%%%%%%%%%%%%%%%
%%%?參數?%%%?
%N?=64;
%fs?=?16;
%M?=?64;
win?=?‘hamming‘;?%?平滑窗類型
if?nargin<1;x=randint（11024）;end;

d_alpha?=?fs/N;?%?1/時間分辨率=循環頻率分辨率
alpha?=?0:d_alpha:30;?%?循環頻率?分辨率=1/時間分辨率
a_len?=?length（alpha）;?%?循環頻率取樣個數

f_len?=?floor（N/M-1）+1;?%?最大平滑窗個數?即頻率采樣個數
f?=?-（fs/2-d_alpha*floor（M/2））?+?d_alpha*M*（0:f_len-1）;?%?頻率采樣點位置

%?S?=?zeros（a_len?f_len）;?%?初始相關功率譜
S?=?zeros（f_lena_len）;?%?初始相關功率譜
i?=?1;?

%%%?信號fft變換?%%%
x?=?fftshift（fft（x（1:N）））;?
x?=x‘;

%%%?遍歷循環頻率取值?%%%
for?alfa?=?alpha

????interval_f_N?=?round（alfa/d_alpha）;?%?循環頻率所對應的頻譜序列序號
????f_N?=?floor（（N-interval_f_N-M）/M）+1;?%?平滑窗的個數

????%%%?生成平滑窗函數?%%%
????g?=?feval（win?M）;?%生成M個點的窗
????window_M?=?g（:?ones（f_N1））;%把g的第一列拿出來存f_n列構成新矩陣，f_n為窗的個數??Mxf_n

????%%%?頻域序列平滑模板?%%%
????t?=?1:M*f_N;
????t?=?reshape（t?M?f_N）;

????%%%?計算X1X2?%%%
????X1?=?x（t）.*window_M;
????X2?=?x（t+interval_f_N）.*window_M;?

????%%%?計算譜相關?%%%
????St?=?conj（X1）.*X2;
????St?=?mean（St?1）;?%?平滑平均
%?????S（i?floor（（f_len-f_N）/2）+（1:f_N））?=?St/N;?%?將結果平移至序列中央以便作圖
????S（floor（（f_len-f_N）/2）+（1:f_N）?i）?=?St/N;?%?將結果平移至序列中央以便作圖
????i?=?i+1;

end

%%%?遍歷循環頻率取值結束?%%%

%%%?循環功率譜作圖?%%%
S=abs（S）./max（max（abs（S）））;
%?mesh（falphaS）;?
figure（4）;
mesh（alphafS）;
%?mesh（S）;
%?axis?tight;
%title（‘BPSK-OFDM‘）;
xlabel（‘循環頻率a‘）;?ylabel（‘頻率f‘）;
title（‘?循環譜‘）

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????1799??2013-11-14?15:32??cs.m