%file name: fsk.m
%Written by Dr. James S. Kang, Cal Poly Pomona.
%Plot of FSK waveform and spectrum.
%d = input data such as [1 1 0 0 1 0 1 0 1 1].
%fb = data rate in bps.
%fm = mark frequency.
%fs = space frequency.
%Ac = Amplitude of the carrier.
%fstart = start frequency for spectrum plot.
%fend = end frequency for spectrum plot.
%To run the program, try
% >>fsk([1 1 0 0 1 0 1 0 1 1],500,1000,2000,1,0,4000)  <Enter>
function[s] = fsk(d,fb,fm,fs,Ac,fstart,fend)
N=size(d,2);	%N = number of data bits.
M=128;	%M = number of samples per bit duration.
tb=1/fb;tm=1/fm;ts=1/fs;
Nm=floor(M*tm/tb); %Nm = number of samples per mark period.
Ns=floor(M*ts/tb); %Ns = number of samples per space period.
step=tb/M;	%step = sampling interval.
for j = 1:N
	if d(j) == 1
	   for i = 1:M
		m((j-1)*M+i)=1;
		s((j-1)*M+i)=Ac*cos(2*pi*(i-1)/Nm);
	   end
	else
	   for i = 1:M
		s((j-1)*M+i)=Ac*cos(2*pi*(i-1)/Ns);
		m((j-1)*M+i)=0;
	   end
	end
end
for k=1:M*N
	t(k)=(k-1)*step;
end
subplot(2,1,1)
plot(t,m)
grid
%xlabel('Time')
ylabel('Amplitude')
title('Input Waveform')
subplot(2,1,2)
plot(t,s)
grid
xlabel('Time')
ylabel('Amplitude')
title('FSK Waveform')
%
%Spectrum of FSK
%
M1=64;	%M1 = number of points per lobe.
fstep=fb/M1; %fstep = freq. step.  
Nf=floor((fend-fstart)/fstep); %Nf = total number of freq. points computed.
for i = 1:Nf
 f(i)=(i-1)*fstep+fstart;
 if f(i) == fm
  S(i)=(Ac^2*tb/8)*(1+(sin(pi*(f(i)-fs)*tb)/(pi*(f(i)-fs)*tb))^2);
 elseif f(i) == fs 
  S(i)=(Ac^2*tb/8)*((sin(pi*(f(i)-fm)*tb)/(pi*(f(i)-fm)*tb))^2+1);
 else
 S(i)=(Ac^2*tb/8)*((sin(pi*(f(i)-fm)*tb)/(pi*(f(i)-fm)*tb))^2+(sin(pi*(f(i)-fs)*tb)/(pi*(f(i)-fs)*tb))^2);
 end
end
S(floor(fm/fstep)+1)=S(floor(fm/fstep)+1)+Ac^2/8;
S(floor(fs/fstep)+1)=S(floor(fs/fstep)+1)+Ac^2/8;
for i = 1:Nf
 dBS(i) = 10*log(S(i));
 if dBS(i) < -200 dBS(i)=-200; end
end
figure
subplot(2,1,1)
plot(f,S)
grid
%xlabel('Frequency')
ylabel('Spectrum')
title('FSK Spectrum (Linear)')
subplot(2,1,2)
plot(f,dBS)
grid
xlabel('Frequency')
ylabel('Spectrum')
title('FSK Spectrum (dB)')
end