% file amfm2.m
% revised from amfmdemo.m
% add m_t to plots
%
% This program illustrates how sidebands
% arise when a carrier wave is modulated.

% The sidebands can be seen on the FFT plots
% and heard by listening to the sound

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

% To select the modulation type, 
% comment out all other equations
% This version is set up for AM

% 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
% (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; 

% select fmod (modulating frequency)
% fixed value or ramp from fmin to fmax
% if want ramp
fmin=8;
fmax=32;
fmod=fmin+(1:fs*T)*fmax/(fs*T);
% if want fixed value
%fmod=16; % suggested values: 16, 64

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

% select beta (modulation index) (=ka for AM)
% fixed value or ramp from betamin to betamax
% if want ramp
betamin=0.5;
betamax=1.5;
%beta=betamin+(1:fs*T)*betamax/(fs*T);
% if want fixed value
beta=1.5; % suggest 0.5, 0.9, 1.0, 1.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;

% EQUATIONS 

% equation for modulating wave m(t)
m_t=cos(twopi_fmod_t);
% general equation for any modulation
% s_t = a_t*cos(twopi_fc_t + phi_t)
%     = x_t*cos(twopi_fc_t) - y_t*sin(twopi_fc_t)

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

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

% 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);
%s_t=x_c+x_s;

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

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

% FFT

% now take fft of s_t
% first divide s_t into blocks of length blockfft
y=reshape(s_t,blockfft,fs*T/blockfft);
ym=reshape(m_t,blockfft,fs*T/blockfft);
% y has blockfft rows, i.e. each column has 
% blockfft elements
nblocks=fs*T/blockfft; % # fft blocks in file
% # cycles of carrier wave per block 
% = fc*T/nblocks = fc*blockfft/fs
% # cycles of modulating wave in file
% = fmod*T/nblocks = fmod*blockfft/fs

z=fft(y); % fft operation applied to each column

% generate plots

figure(1); % 2-D view of FFTs
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); % 3-D view of FFTs
% 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));hold on
plot((1:blockfft)/fs,ym(:,1),'r-');hold off;
title('first block (time domain)');
xlabel('time (sec)');
ylabel('amplitude');

figure(4); %time domain plot of middle block
plot((1:blockfft)/fs,y(:,nblocks/2));hold on
plot((1:blockfft)/fs,ym(:,nblocks/2),'r-');hold off;
title('middle block (time domain)');
xlabel('time (sec)');
ylabel('amplitude');

figure(5); %time domain plot of last block
plot((1:blockfft)/fs,y(:,nblocks));hold on
plot((1:blockfft)/fs,ym(:,nblocks),'r-');hold off;
title('last block (time domain)');
xlabel('time (sec)');
ylabel('amplitude');

% output sound and .wav file

s_tmax=max(s_t);
soundsc(s_t/s_tmax,fs);
wavwrite(s_t/s_tmax,fs,'amfm2.wav');
%auwrite(s_t/s_tmax,fs,'amfm2.wav'); for unix