T=1;

n=1;
m=2;


f=0:0.03:5.1;
N=171;
y= zeros(1,N);
z=zeros(1,N);
phi0= 16*T*T/pi/pi;
for j=1:N
phi_msk = 16*T/(pi*pi)*(cos(2*pi*f(j)*T)^2 )/((16*T*T*f(j)*f(j)-1)^2) ;
phi_fsk = T/4/pi/pi* ( ( (n*sin(2*pi*f(j)*T-n*pi)/(4*f(j)*f(j)*T*T - n*n)) - (m*sin(2*pi*f(j)*T-m*pi)/(4*f(j)*f(j)*T*T- m*m))  )^2 ) ;
y(j)= phi_msk;
z(j)= phi_fsk;
end

plot (f,y, f,z,'--' );
legend('MSK','BFSK');