%%%  設(shè)置調(diào)制參數(shù)  bpsk qpsk msk
fc=200;
fd=40;
fs=1000;

%%%  生產(chǎn)帶通信號(hào)

% s=randint(1,400,[0,1]);
% y=dmod(s,fc,fd,fs,'psk',2);
% h1=fir1(50,[(fc-fd)/fs*2,(fc+fd)/fs*2]);
% y=conv(y,h1);

% s=randint(1,400,[0,3]);
% y=dmod(s,fc,fd,fs,'psk',4);
% h1=fir1(50,[(fc-fd)/fs*2,(fc+fd)/fs*2]);
% y=conv(y,h1);

% s=randint(1,400,[0,1]);
% y=dmod(s,fc,fd,fs,'msk');

%%% 設(shè)置調(diào)制參數(shù) am fm
fc=200;
fs=600;
%%%  生產(chǎn)帶通信號(hào)
n=randn(1,2200);
Xm(1)=n(1);
for i=2:2200
    Xm(i)=0.979274*Xm(i-1)+n(i);
end
y=amod(Xm,fc,fs,'amdsb-sc');
% y=amod(Xm,fc,fs,'fm');

%%%  生產(chǎn)解析信號(hào)
h=hilbert(y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  產(chǎn)生用于dsp調(diào)試的數(shù)據(jù)
N=length(h);
fid=fopen('input_data.dat','w');
fprintf(fid,'1651 4 0 0 0\r\n');
for i=1:N
    fprintf(fid,'%f\r\n',real(h(i)));
    fprintf(fid,'%f\r\n',imag(h(i)));
end
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  計(jì)算循環(huán)自相關(guān)  2048點(diǎn)fft
spec=[];
J=zeros(1,2048);
M=10;
x=h(129:129+2047);
for tao=0:M-1
  xn=h(129-tao:129-tao+2047);
  z=x.*conj(xn);           % 共軛相乘
  spec1=abs(fft(z)).^2;    % 計(jì)算頻譜
  J=J+spec1;               % 頻譜累加
  spec=[spec;spec1];       % 存儲(chǔ)頻譜，用于調(diào)試時(shí)觀測(cè)
end

%%% 畫循環(huán)自相關(guān)
subplot(2,1,1);
mesh(spec(:,2:1024));      % 為了觀測(cè)方便，去掉第一列（對(duì)應(yīng)的循環(huán)頻率為0）
subplot(2,1,2);
plot(J(2:1024));           % 為了觀測(cè)方便，去掉第一列（對(duì)應(yīng)的循環(huán)頻率為0）

%%% 識(shí)別
thro1=1;                   % 判決門限
JJ=J(1:1024);JJ(1)=0;      % 找到2--1024之間最高的譜線
[val alpha]=max(JJ) ;      
                           

% %%% 第一種判決規(guī)則
% N=5;                       % 取最高譜線左右兩邊各N根譜線之和做歸一化
% t=[alpha,alpha+1:alpha+N];       % 找到用于歸一化的下標(biāo)
% for i=1:N
%     if alpha-i==0
%         temp=2048;
%     else
%         temp=mod(alpha-i,2048);
%     end
%     t=[temp t];
% end
% num=sum(J(t(5:7)));
% den=sum(J(t))-num;
% Rx=num/den                 % 歸一化
% esti_fd=(alpha-1)*fs/2048;
% 
% if Rx < thro1
%     str='it is am,fm or msk' 
% else
%     str='it is bpsk or qpsk'
%     esti_fd=(alpha-1)*fs/2048
% end


%%% 第二種判決規(guī)則
if alpha<10+1+1
%     str='it is am,fm or msk'
    env=abs(h(1:2048)).^2;
    m_env=mean(env);
    env_cn=(env/m_env)-1;
    env_fft=abs(fft(env_cn)).^2;    % 此處不應(yīng)該平方，應(yīng)該除以2048，為了降低dsp計(jì)算復(fù)雜度，省略此除法運(yùn)算
    env_fft(1)=0;
    max_env=max(env_fft);
    if max_env<1000
        str='it is fm or msk'
    else
        str='it is am'
    end
else
    t=[alpha-10:alpha+10];  % 生成用于歸一化的下標(biāo)
    num=sum(J(t(10:12)));
    den=sum(J(t))-num;
    Rx=num/den;                 % 歸一化
    if Rx < thro1
%         str='it is am,fm or msk' 
        env=abs(h(1:2048)).^2;
        m_env=mean(env);
        env_cn=(env/m_env)-1;
        env_fft=abs(fft(env_cn)).^2;    % 此處不應(yīng)該平方，應(yīng)該除以2048，為了降低dsp計(jì)算復(fù)雜度，省略此除法運(yùn)算
        env_fft(1)=0;
        max_env=max(env_fft);
        if max_env<1000
            str='it is fm or msk'
        else
            str='it is am'
        end
    else
        esti_fd=(alpha-1)*fs/2048
        h2fft=abs(fft(h(1:2048).^2));         % 平方，F(xiàn)FT
        m_h2fft=mean(h2fft);
        h2fft_cn=(h2fft/m_h2fft)-1;
        [max_h2fft index]=max(h2fft);
        num=sum(h2fft(index-2:index+2));
        den=sum(h2fft(index-10:index+10))-num;
        rfc=num/den;
        if rfc>1
            esti_fc=(index-1)/2048*fs/2
            str='it is bpsk'
        else
            h4fft=abs(fft(h(1:2048).^4));         % 平方，F(xiàn)FT
            [max_h4fft index]=max(h4fft);
            esti_fc=(index-1)/2048*fs/4
            str='it is qpsk'
        end
    end
end



h2fft=abs(fft(h(1:2048)));         % 平方，F(xiàn)FT
plot(h2fft)
% m_h2fft=mean(h2fft);
% h2fft_cn=(h2fft/m_h2fft)-1;
% [max_h2fft index]=max(h2fft);
% num=sum(h2fft(index-2:index+2));
% den=sum(h2fft(index-10:index+10))-num;
% rfc=num/den
% 
% % var(h2fft_cn)
% plot(h2fft_cn)