% Plot of FM power spectra
%
clf
DELfTb_min = input('Enter min. freq spacing times bit period between tones ');
DELfTb_0 = input('Enter step size in frequency spacing X bit period ');
fTb = 0:.01:10;
fcTb = 5;
for n=1:6
   DELfTb = DELfTb_min + (n-1)*DELfTb_0
   f1Tb = fcTb-DELfTb/2;
   f2Tb = fcTb+DELfTb/2;
   alpha = 0.5*(f1Tb + f2Tb);
   beta = 0.5*(f2Tb - f1Tb);
   num_plus = ((sin(pi*(fTb+f1Tb+eps*pi))).^2).*(sin(pi*(fTb+f2Tb+eps*pi))).^2;
   num_minus = ((sin(pi*(fTb-f1Tb+eps*pi))).^2).*(sin(pi*(fTb-f2Tb+eps*pi))).^2;
   den_plus = 2*pi^2*(1-2*cos(2*pi*(fTb+alpha+eps*pi)).*cos(2*pi*beta+eps*pi)+(cos(2*pi*beta+eps*pi)).^2);
   den_minus = 2*pi^2*(1-2*cos(2*pi*(fTb-alpha+eps*pi)).*cos(2*pi*beta+eps*pi)+(cos(2*pi*beta+eps*pi)).^2);
   term_plus = (1./(fTb+f1Tb+eps*pi) - 1./(fTb+f2Tb+eps*pi)).^2;
   term_minus = (1./(fTb-f1Tb+eps*pi) - 1./(fTb-f2Tb+eps*pi)).^2;
   G_plus = num_plus./den_plus.*term_plus;
   G_minus = num_minus./den_minus.*term_minus;
   G = G_plus+G_minus;
   area = sum(G)*.01            % Check on area under spectrum
   GN = G/area;
   subplot(3,2,n),xlabel('fT_b'),plot(fTb, GN),ylabel('Spectral level'),axis([0 10 0 max(GN)]),...
   legend(['\deltafT_b=' ,num2str(DELfTb)]),...
end