%%%% 本程序思想：初步認定為頻域加窗的FFT（FAM）算法，運用信號與其復共扼頻譜相乘的方式來表達循環譜。處理信號為運用BPSK調制產生的
%數據設置部分
Rb=2;     %碼源速率
Ts=1/Rb;   %碼源間隔
L=2^6;     %單位時間內信號抽樣點數
num=16;   %一個段落中碼源個數
N=num*L;  %一個段落內總的抽樣點數 (1024)

f0=16;   %載波頻率
dt=Ts/L;  %抽樣時間間隔
df=1/(dt*N);   %頻率間隔
Bs=df*N/2;   %帶寬
fs=Bs*2;
f=df/2-Bs:df:Bs;   %抽樣頻率
t=dt/2:dt:num*Ts;  %抽樣時間  t= 1/(4*64) : 1/(2*64) : 8;

a=sign(randn(1,num)); %碼源序列 generate the initial signal (length is 16)

%%%%%%%%%%%% BPSK調制信號產生
%通過基帶成型濾波器后的信號   function ?????????????????????????????????
for rr=1:num
    I((rr-1)*L+1:rr*L)=a(rr);  % number of I is 1024 .extend the length of signal 'a' to 64 times length.
end

%最終調制信號 the signal modulated by BPSK carrier
s=I.*cos(2*pi*f0*t);  % the number of 't'and 's' are the same.(1024)
% plot(s,'ro');
% axis([0 35 -1 1]); look 35/1024 percent of the plot.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% begin to process the signal (look the .doc )
M=2; %  ??????????????????????????
Fs=fs;
N = (M*Fs)/df;  % the new value of 'N'is 2 times the value of old 'N'. (2048)
N = pow2 (nextpow2(N)); % windowing record for FFT (2048)   % NEXTPOW2(N) returns the first P such that 2^P >= abs(N)
                                                      % X = POW2(Y) for each element of Y is 2 raised to the power Y
%設置FFT點數為1024
X = fft(s,N);   % fft of the truncated(被刪節的) (or zero padded) time series (2048 columns)
X = fftshift(X);% shift components of fft  (2048 columns)
Xc = conj(X);                 % precompute the complex conjugate(復共軛) vector (2048 columns)

S = zeros (N,N);              % size of the Spectral Correlation Density matrix
f = zeros (N,N);              % size of the frequency matrix;
alfa = zeros (N,N);           % size of the cycle frequency matrix
F = Fs/(2*N);                 % precompute(預先計算的) constants -  F = Fs/(2*N);     
G = Fs/N;                     % precompute constants -  G = Fs/N;  
m = -M/2+1:M/2;               % set frequency smoothing window index  % m(0)=0,m(1)=1.

for k = 1:N                                % fix k    (2048 columns)
    % computes vectors of f and alfa,
    % store frequency and cycle frequency data for given k.
    k1 = 1:N;
    f(k,k1) = F*(k+k1-1) - Fs/2;          % Computes f values and shift them to center in zero (f = (K+L)/2N) [1]
    alfa(k,k1) = G*(k-k1 + N-1) - Fs;       % Computes alfa values and shift them to center in zero (alfa = (K-L)/N) [1]
       for k1 = 1:N                          %fix k1 = J
            %calculate X(K+m) & conj (X(J+m)) for arguments of X(1:N) only
           B = max(1-k, 1-k1);          % Largest min of 1 <= (K+m)| (J+m) <= N
           A = min (N-k, N-k1);         % Smallest max of 1 <= (K+m)| (J+m) <= N
           n = m((m<=A) & (m>=B));      % fix the index out of range problem by truncating the window.  % m(0)=0,m(1)=1.
             if isempty(n)
                S(k,k1) = 0;
             else
                p = k+n; q = k1+n;
                Y = X(p).*Xc(q);  
                S(k,k1) = sum(Y);  %  the power spetral Correlation destiny function ??????????????????
             end
       end
end
% haha=length(sign([real(Y)]))
S = abs(S./max(max(S)));% normalize output matrix
% figure(1);
% contour (alfa, f, S); grid;
figure(1);
f=zeros(N,N);
plot(alfa,S,'r'); % 'alfa'is cycle frequency.'f'is frequency.'S'is Spetral Correlation Density Function.
axis([0 40 0 1]);
ylabel('Frequency (Hz)');xlabel('Cycle frequency (Hz)'); 
title('BPSK循環譜三維圖');