% 
% 2fsk.m
% Simulation program to build 2FSK
% All of the programs had been built successfully in MATLAB 6.5

%************************ Preparation part *****************************
NUM = 1024;            % The number of FFT point
sr=1000000.0;          % Symbol rate
fc = 1000000;          % The frequency of the carrier
fs = 32000000;         % The sample rate
ml=1;                  % Number of modulation levels
br=sr.*ml;             % Bit rate (=symbol rate in this case)
nd=NUM*fc/fs;         % The number of symbols
n =fc/sr;              % The number of the carrier in a symbol
ebno=10;               % Eb/N0
IPOINT=32;             % Number of oversamples
Ts = 1/sr;             % The period of the sourse signal

t=0:1/fs:1/fc*nd*n-1/fs;    % The sample point in one carrier period
fcos = cos(2*pi*fc*t);      % cosin carrier
kf=800000;
data5(1)=0;

%*********************** Filter initialization **************************

irfn=21;                                 % Number of filter taps
alfs=0.33;                                % Rolloff factor
[xh]=hrollfcoef(irfn,IPOINT,sr,alfs,1);  % Transmitter filter coefficients

%*********************** Data generation **********************************

data=rand(1,nd)>0.5;                     % rand: built in function

%*********************** 2fsk modulation **********************************

data1=data.*2-1;
[data2]=oversample(data,nd,IPOINT);
data3=conv(data2,xh);                                %(conv: build in function)
data4=data3(irfn*IPOINT/2-1:NUM+irfn*IPOINT/2-2);    % date4 is the baseband singal of 2fsk 
data4=data4+0.3;                                     
data5(1)=data4(1)*(1/fs);
for i=1:length(data4)-1;                              
    data5(i+1)=data5(i)+data4(i)*(1/fs);
end

modulation = cos(2*pi*fc*t+2*pi*kf*data5);                % The modulation is the 2fsk singal

%**********************comppute the amplitude *********************************************

mod_hilb=hilbert(modulation);
x=real(mod_hilb);
y=imag(mod_hilb);
z=sqrt(modulation.*modulation+y.*y) ;                  % z is the amplitude

x(1)=0;
for i=1:length(z)-1                             % 
    x(i+1)=x(i)+z(i);
     ma=x(1024)/NUM;
end

for i=1:length(z)-1 
  x(i)=z(i)./ma;
  x(i)=x(i)-1;
  x(i)=x(i)*x(i);
end
  y=fft(x,NUM);
  yxk=abs(y(1:NUM/2));
 for i=1:NUM/2
    for j=1:NUM/2-i
        t=yxk(j);yxk(j)=yxk(j+1);yxk(j+1)=t;
    end
end
max=yxk(NUM/2)/NUM;
fprintf('\n%f',max);