% file amfmdemo.m

%(To save this file to run on Matlab,
% use Save As in the Browser, and save as
% amfmdemo.m)

%
% This program illustrates how sidebands
% arise when a carrier wave is modulated.

% This version includes equations for
% all modulation schemes
% AM, DSB  (ring), SSB, FM, PM

% To select the modulation type, 
% comment out all other equations.
% In this version, all but the AM equation
% is commented out.

% Consider 256 Hz carrier wave
% sampled at 8192 Hz, i.e. about 32 samples
% per cycle.
% Modulate this carrier with 
% modulation frequency fmod, and
% modulation index beta.
% Vary either fmod or beta linearly
% from near 0 to a maximum value.
% Observe spectrum and listen to sound.
% Select parameters by editing this file.
% Future versions may have a menu and/or GUI,
% and may have real time spectrum plots
% in sync with the audio.

fs=8192; % sampling rate, suggest 8192 (power of 2)
T=16; %wavefile time in sec, suggest 16
fc=256; % carrier frequency, suggest 256 (middle C)
blockfft=1024; % fft block size, suggest 1024
% frequency resolution f0=fs/N0 = fs/blockfft
% suggest f0 = 8192/1024 = 8 Hz

% want to write cos (2 pi fc t)
% but in sampled systems t = n*ts = n/fs
% where n is an integer index and fs=1/ts is the 
% sampling rate 
% fs*T is the total number of samples in the file
% # samples = cycles*fs/f0 = f0*T*fs/f0 = fs*T

% (1:fs*T)/(fs*T) is vector of fs*T numbers
% ranging from 0+ to 1
% index n is represented by vector of increasing
% integers (1:fs*T) = [ 1 2 3 ... fs*T]

% 2 pi fc t is written as below
twopi_fc_t=(1:fs*T)*2*pi*fc/fs; 

% vary fmod
fmin=1;
fmax=16;
%fmod=fmin+(1:fs*T)*fmax/(fs*T);
fmod=16; % suggested values: 16, 64

% 2 pi fmod t is written as below
twopi_fmod_t=(1:fs*T)*2*pi.*fmod/fs; 

% vary beta (modulation index)
betamin=0;
betamax=2;
beta=betamin+(1:fs*T)*betamax/(fs*T);
%beta=0.5; % suggest 0.5, 1, 2, 4, 8

% for FM only
% vary delf (frequency deviation)
% beta = delf/fmod or delf=beta*fmod
delfmin=0;
delfmax=64;
delf=delfmin+(1:fs*T)*delfmax/(fs*T);
%delf=32; % suggest 8, 16, 32, 64
% for FM only, not for PM
% beta=delf/fmod;

% equation for modulating wave m(t)
m_t=cos(twopi_fmod_t);

% equation for AM signal x(t)
x_t=(1+beta.*cos(twopi_fmod_t)).*cos(twopi_fc_t);

% equation for DSB (ring modulation)
%x_t=(beta.*cos(twopi_fmod_t)).*cos(twopi_fc_t);

% equation for SSB 
%x_c=(beta.*cos(twopi_fmod_t)).*cos(twopi_fc_t);  
%x_s=(beta.*sin(twopi_fmod_t)).*sin(twopi_fc_t);
%x_t=x_c+x_s;

% equation for FM
%x_t=cos((twopi_fc_t)+beta.*sin(twopi_fmod_t));

% equation for PM
%x_t=cos((twopi_fc_t)+beta.*cos(twopi_fmod_t));

% now take fft of x_t
% first divide x_t into blocks of length blockfft
y=reshape(x_t,blockfft,fs*T/blockfft);
% y has blockfft rows, i.e. each column has 
% blockfft elements
nblocks=fs*T/blockfft; % # fft blocks in file
z=fft(y); % fft operation applied to each column

figure(1);
za=z(:);maxz=max(abs(za));% find max of all ffts
%plot((1:blockfft)*fs*T/blockfft,abs(z(:,1)));
plot((1:blockfft)*fs/blockfft,abs(z));
axis([fc-2*fmax fc+2*fmax 0 maxz]);
title('2-D view of FFTs');
xlabel('frequency (Hz)');
ylabel('amplitude');

figure(2);
% set up xx and yy axes for 3-d mesh plot
xx=(1:blockfft)*fs/blockfft;
yy=(1:fs*T/blockfft);
zt=z';
mesh(xx,yy,abs(zt)); % makes 3-d plot
axis([fc-2*fmax fc+2*fmax 0 fs*T/blockfft 0 maxz]);
view(25,30);
title('3-D view of FFTs');
xlabel('frequency (Hz)');
zlabel('amplitude');
ylabel('block #');

figure(3); %time domain plot of first block
plot((1:blockfft)/fs,y(:,1));
title('first block (time domain)');
xlabel('time (sec)');
ylabel('amplitude');


figure(4); %time domain plot of last block
plot((1:blockfft)/fs,y(:,nblocks));
title('last block (time domain)');
xlabel('time (sec)');
ylabel('amplitude');


x_tmax=max(x_t);
soundsc(x_t/x_tmax,fs);
wavwrite(x_t/x_tmax,fs,'amfmdemo.wav');
%auwrite(xout/10,48000,'clar13.au');