%% 編一個 4fsk 的調制方式的程序；
%% 基帶信號是一串隨機碼，每個基帶碼元用一個周期的載波發送；
%% 設置采樣頻率，每個碼元的采樣點數，載波頻率；
function [sig_M]=fsk_4(M)
         %% 每一個碼元發送一個周期的載波
         fs=1000; N=50; f=20; t=(1:N)'/fs; 
%          sig=sin(2*pi*f*t);
         % 設置發送的原始隨機碼序列
         M=M; Ji=rand(M,1); Ji=ceil(Ji-0.5); 
         %% 把原來的2進制序列轉換為4進制序列,長度會減半；
         MM=M/2; Ji_2=zeros(MM,1);
         for i=1:MM
             if (Ji(2*i-1)==0)&&(Ji(2*i)==0)
                 Ji_2(i)=0;
             elseif (Ji(2*i-1)==0)&&(Ji(2*i)==1)
                 Ji_2(i)=1;
             elseif (Ji(2*i-1)==1)&&(Ji(2*i)==0)
                 Ji_2(i)=2;
             elseif (Ji(2*i-1)==1)&&(Ji(2*i)==1)
                 Ji_2(i)=3;
             end;
         end;        
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %% 設計調制信號，共有4種頻率,每個碼元的發送時間和2ask是一樣的；
   %% sig_M表示調制后的信號;
      sig_M=zeros(MM*N,1);
for i=1:MM
    if Ji_2(i)==0
        sig_M((i-1)*N+1:1:(i*N))=sin(2*pi*f*t);
    elseif Ji_2(i)==1
        sig_M((i-1)*N+1:1:(i*N))=sin(2*pi*(f+f)*t);
    elseif Ji_2(i)==2
        sig_M((i-1)*N+1:1:(i*N))=sin(2*pi*(f+2*f)*t);
    elseif Ji_2(i)==3
        sig_M((i-1)*N+1:1:(i*N))=sin(2*pi*(f+3*f)*t);
    end;
end;
     sig_M=sig_M;
     %%%% 開始畫圖；
     figure;
     subplot(321);
     stem(Ji);title('2進制基帶時域波形');grid;axis tight;

     subplot(322);
     stem(Ji_2);title('4進制基帶時域波形');grid;axis tight;

     subplot(323);
     plot(sig_M);title('4fsk時域波形');grid;axis tight;

     sig_fft1=abs(fftshift(fft(sig_M)));
     subplot(324);
     plot(fs*((0:MM*N-1)-(MM*N)/2)/(MM*N),(sig_fft1));
     title('4fsk頻譜');grid;axis tight;

     sig_noise=awgn(sig_M,10,'measured');
     subplot(325);
     plot(sig_noise);title('4fsk時域噪聲波形');grid;axis tight;

     sig_fft2=abs(fftshift(fft(sig_noise)));
     subplot(326);
     plot(fs*((0:MM*N-1)-(MM*N)/2)/(MM*N),(sig_fft2));
     title('4fsk含噪聲的頻譜');grid;axis tight;