% file amfmdem4.m
%
% from amfmdem3.m
% should run on Matlab 4.2, function calls deleted
% modified from file amfm2.m
% by adding receivers and filters
%
% This program illustrates how sidebands
% arise when a carrier wave is modulated.
% It includes both modulators (transmitters)
% and demodulators (receivers) for 
% modulation types AM, DSB (ring), SSB, FM, PM, 
% combined AM/FM and general i-q (AM/PM).

% The program includes various low pass, bandpass
% and notch filters.
% Output includes the sound of modulated signals,
% time and frequency domain plots of modulated
% and demodulated signals.

% To select the modulation type, and receiver 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 betaf.
% Vary either fmod or betaf 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.

% =======================================
% MATLAB HELP
% to learn more about Matlab, type 'intro'
% (without quotes) at the command prompt
% also type 'demo' for more complete intro
% for help on any function, type 'help functionname'
% for help on syntax:
% type 'help punct' 'help ops'  'help arith'
% for web-based help, type 'helpdesk'
% keyword search, type 'lookfor keyword'
% to retain variables in workspace on error
% type 'dbstop if error' before running 

% made amfmdem3.m file into a function
% so can have other functions
% to keep values of variables in workspace
% for viewing after running
% need to activate keyboard function
% at end of main routine
% before the other functions
% then when finished running and viewing
% get K at the command line promt
% type 'return' and then 'clear all'
% before running again
% This version includes equations for
% all modulation schemes
% AM, DSB (ring), SSB, FM, PM


% =======================================
% set universal parameters
% =======================================
fs=8192; % sampling rate, suggest 8192 (power of 2)
T=2; %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

% =======================================
% TRANSMITTER (generate modulated signal) 
% =======================================

% =======================================
% preliminaries: - write cosine waves
% and set parameters fmod, betaf, delf
% =======================================

% 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]
% for more info, type 'help colon'

% carrier wave
% 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-fmin)/(fs*T);
% if want fixed value
%fmod=32; % suggested values: 16, 64
% 2 pi fmod t is written as below
twopi_fmod_t=(1:fs*T)*2*pi.*fmod/fs; 

% select betaf (modulation index) (=ka for AM)
% fixed value or ramp from betamin to betamax
% if want ramp
betamin=0;betamax=2;
%betaf=betamin+(1:fs*T)*(betamax-betamin)/(fs*T);
% if want fixed value
betaf=1; % suggest 0.5, 0.9, 1.0, 1.1, 2, 4, 8

% for FM only
% vary delf (frequency deviation)
% set only one of delf, betaf, then
% betaf = delf/fmod or delf=betaf*fmod
delfmin=8;
delfmax=32;
%delf=delfmin+(1:fs*T)*(delfmax-delfmin)/(fs*T);
delf=16; % suggest 8, 16, 32, 64
% for FM only, not for PM, if set delf, then
% betaf=delf/fmod;

% =========================
% MODULATION EQUATIONS
% =========================
% for each modulation type, there may be
% one or more equivalent equations

% equation for modulating wave m(t)
m_t=cos(twopi_fmod_t);
% equation for carrier wave c(t)
c_t=cos(twopi_fc_t);

% general equation for any modulation
% first choose mapping from m(t) to 
% (a_t,phi_t) or (x_t,y_t) for the 
% desired modulation
% s_t = a_t*cos(twopi_fc_t + phi_t)
% s_t = x_t*cos(twopi_fc_t) - y_t*sin(twopi_fc_t)

% equation for AM signal s(t)
s_t=(1+betaf.*cos(twopi_fmod_t)).*cos(twopi_fc_t);
% a_t=1+betaf.*cos(twopi_fmod_t);
% phi_t=0;
% s_t=a_t.*cos(twopi_fc_t);
% s_t=(1+betaf.*m_t).*c_t;

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

% equation for SSB 
%x_c=(betaf.*cos(twopi_fmod_t)).*cos(twopi_fc_t);  
%x_s=(betaf.*sin(twopi_fmod_t)).*sin(twopi_fc_t);
%s_t= x_c - x_s; %USB
% alternate method for SSB
% using hilbert transform mhat_t=HT(m_t)
%mhat_t = imag(hilbert(m_t)); % sin(twopi_fmod_t)
%cq_t = imag(hilbert(c_t)); % sin(twopi_fc_t);
%s_t=m_t.*c_t -mhat_t.*cq_t; %USB

% equation for FM
%s_t=cos((twopi_fc_t)+betaf.*sin(twopi_fmod_t));
%s_t=cos((twopi_fc_t)+(delf./fmod).*sin(twopi_fmod_t));
%s_t=cos((twopi_fc_t)+phi_t);
% where
% phi_t = 2*pi*delf * integral(m_t)dt
% integral(m_t) = cumsum(cos(twopi_fmod_t)/fs)
%       = sin(twopi_fmod_t)/(2*pi*fmod)
% assuming no DC term in m_t
% s_t=cos((twopi_fc_t)+2*pi*delf.*cumsum(m_t)/fs);

% equation for combined AM/FM
% with arbitrary functions for 
% amplitude envelope a_t, and
% modulation index envelope b_t (time-varying betaf)
% often used for FM synthesis of instrument sounds
% for bell sound, try fmod = 2*fc, b_t(1)=10
%tau1=T/3; % decay constant relative to file length
%tau2=T/2;
% a_t=exp(-(1:fs*T)/(fs*tau1));
% b_t=10*exp(-(1:fs*T)/(fs*tau2)); % b_t(1)=10
% s_t=a_t.*cos((twopi_fc_t)+b_t.*sin(twopi_fmod_t));

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

% impulse signal
%s_t=[1 zeros(1,fs*T-1)];

% equations for addition of two cos waves
% at fc+fmod and fc-fmod
% the sum will be the same as a DSB signal
% with carrier fc and modulation fmod 
%fcplusfmod=fc+fmod;
%fcminusfmod=fc-fmod;
% 2 pi fcplusfm t is written as below
%twopi_fcplusfmod_t=(1:fs*T)*2*pi.*fcplusfmod/fs; 
% 2 pi fcminusfm t is written as below
%twopi_fcminusfmod_t=(1:fs*T)*2*pi.*fcminusfmod/fs; 
% equation for the sum and difference waves
%csum_t=cos(twopi_fcplusfmod_t);
%cdif_t=cos(twopi_fcminusfmod_t);
% equation for beats between two waves
% s_t=csum_t+cdif_t;
% s_t is equivalent to DSB with fc and fmod

% =======================================
% FILTERS
% =======================================
% Set filter coefficient vectors for
% 3 filters, LPF fmax, BPF fc and 2fc.
% There are two vectors for each filter:
% b (numerator) and a (denominator) of H(z)
% then use Matlab filter function
% with input vector x to get output
% y = filter(b,a,x) 
% =======================================
% filter gains and bandwidths may need
% adjustment

% low pass filter, 1 pole 1 zero, cutoff fmax
% gain set to be unity at zero frequency (z=1)
% H(z)=gain*(z+1)/(z-p)
% pole p on real axis, zero at -1
% magnitude pmag=|p| near but < 1
% choose pole pmag for about 3dB down at fmax
pmag=1-(1/2)*(2*pi*fmax/fs);
gain=(1-pmag)/2; % to get unity gain at z=1
bfmax=gain*[1 1];% numerator (zeros)
afmax=[1 -pmag];% denominator (poles)
%plotLPFfmax(afmax,bfmax,fmax,fc,fs);
%function plotLPFfmax(a,b,fmax,fc,fs);
% plot frequency response
a=afmax;
b=bfmax;
zz=exp(j*(0:fs/2-1)*2*pi/fs); % z values along unit circle
H=(b(1)*zz+b(2))./(a(1)*zz+a(2));
%H(H==0)=10^(-6); % avoid H=0 to avoid log of zero
H(1)=10^(-6);
figure(9);
subplot(2,1,1);
plot((1:fs/2)-1,abs(H));
axis([0 2*fmax 0 1]);
subplot(2,1,2);
plot((1:fs/2),20*log10(abs(H/max(H))));
axis([0 2*fc -40 0]);
title('|H(z)| 1-pole, 1-zero LPF');
xlabel('frequency (Hz)');
ylabel('amplitude');

% low pass filter, 2 pole 2 zero, cutoff fmax
% gain set to be unity at zero frequency (z=1)
% H(z)=gain*(z+1)(z+1)/(z-p)(z-p)
%     =gain*(z^2+2*z+1)/(z^2-2*p*z+p^2)
% double pole p on real axis, double zero at -1
% magnitude pmag=|p| near but < 1
% choose pole pmag for about 3dB down at fmax
% 2-sided 3dB BW(Hz) = (1-pmag)/(pi*sqrt(pmag))
pmag2=1-(2*pi*fmax/fs);
gain2=(1-2*pmag2+pmag2^2)/4; % to get unity gain at z=1
bfmax2=gain2*[1 2 1];% numerator (zeros)
afmax2=[1 -2*pmag2 pmag2^2];% denominator (poles)
%plotLPFfmax2(afmax2,bfmax2,fmax,fc,fs);
%function plotLPFfmax2(a,b,fmax,fc,fs);
% plot frequency response
a=afmax2;
b=bfmax2;
zz=exp(j*(0:fs/2-1)*2*pi/fs); % z values along unit circle
H=(b(1)*zz.^2+b(2)*zz+b(3))./(a(1)*zz.^2+a(2)*zz+a(3));
%H(H==0)=10^(-6); % avoid H=0 to avoid log of zero
H(1)=10^(-6);
figure(7);
subplot(2,1,1);
plot((1:fs/2)-1,abs(H));
axis([0 2*fmax 0 1]);
subplot(2,1,2);
plot((1:fs/2),20*log10(abs(H/max(H))));
axis([0 2*fc -40 0]);
title('|H(z)| 2-pole, 2-zero LPF');
xlabel('frequency (Hz)');
ylabel('amplitude');


% 2-pole 2-zero bandpass filter at fc
% zeros at +1, -1, poles at p=|p|e^+-j(theta)
% where theta = 2pi fc/fs
% for peak response at fc
% H(z) = (z+1)(z-1)/(z-p)(z-p*)
%      = (z^2-1)/(z^2-2*|p|*cos(theta)*z+|p|^2)
theta3=2*pi*(fc)/fs;
p3=0.95; % sets bandwidth
% 2-sided 3dB BW(Hz) = (1-p1)/(pi*sqrt(p1))
%gain= % should scale for unity output at peak
afc=[1 -2*p3*cos(theta3) p3*p3]; % denominator (poles)
bfc=[1 0 -1]; % numerator (zeros)
%plotBPFfc(afc,bfc,fc,fs);
%function plotBPFfc(a,b,fc,fs);
% plot frequency response
a=afc;
b=bfc;
zz=exp(j*(0:fs/2-1)*2*pi/fs); % z values along unit circle
H=(b(1)*zz.^2+b(2)*zz+b(3))./(a(1)*zz.^2+a(2)*zz+a(3));
%H(H==0)=10^(-6); % avoid H=0 to avoid log of zero
H(1)=10^(-6);
figure(10);
subplot(2,1,1);
plot((1:fs/2)-1,abs(H));
axis([0 2*fc 0 max(abs(H))]);
subplot(2,1,2);
plot((1:fs/2),20*log10(abs(H/max(H))));
axis([0 fs/2 -40 0]);
title('|H(z)| 2-pole, 2-zero BPF');
xlabel('frequency (Hz)');
ylabel('amplitude');

% 2-pole 2-zero bandpass filter at 2*fc
% zeros at +1, -1, poles at p=|p|e^+-j(theta)
% where theta = 2pi (2fc)/fs
% for peak response at fc
% H(z) = (z+1)(z-1)/(z-p)(z-p*)
%      = (z^2-1)/(z^2-2*|p|*cos(theta)*z+|p|^2)
theta4=2*pi*(2*fc)/fs;
p4=0.95; % sets bandwidth
% 2-sided 3dB BW(Hz) = (1-p)/(pi*sqrt(p))
%gain= % should scale for unity output at peak
a2fc=[1 -2*p4*cos(theta4) p4*p4]; % denominator (poles)
b2fc=[1 0 -1]; % numerator (zeros)
%plotBPF2fc(a2fc,b2fc,fc,fs);
%function plotBPF2fc(a,b,fc,fs);
% plot frequency response
a=a2fc;
b=b2fc;
zz=exp(j*(0:fs/2-1)*2*pi/fs); % z values along unit circle
H=(b(1)*zz.^2+b(2)*zz+b(3))./(a(1)*zz.^2+a(2)*zz+a(3));
%H(H==0)=10^(-6); % H(1)=0 is really zero, want to avoid log(0)
H(1)=10^(-6);
figure(8);
subplot(2,1,1);
plot((1:fs/2)-1,abs(H));
axis([0 4*fc 0 max(abs(H))]);
subplot(2,1,2);
plot((1:fs/2),20*log10(abs(H/max(H))));
axis([0 fs/2 -40 0]);
title('|H(z)| 2-pole, 2-zero BPF');
xlabel('frequency (Hz)');
ylabel('amplitude');

% 2-pole 2-zero notch filter at fc
% zeros at z1=e^+-j(theta),
% poles at p=|p|e^+-j(theta)
% where theta =  2*pi fc/fs
% for peak response at fc
% H(z) = (z-z1)(z-z1*)/(z-p)(z-p*)
%      = (z^2-(cos(theta)+1)/(z^2-2*|p|*cos(theta)+|p|^2)
theta5=2*pi*(fc)/fs;
p5=0.95; % sets bandwidth
% 2-sided 3dB BW(Hz) = (1-p1)/(pi*sqrt(p1))
gain5= (1-2*p5*cos(theta5)+p5^2)/(2-2*cos(theta5)); 
anfc=[1 -2*p5*cos(theta5) p5*p5]; % denominator (poles)
bnfc=gain5*[1 -2*cos(theta5) 1]; % numerator (zeros)
%plotBPFnfc(anfc,bnfc,fc,fs);
%function plotBPFnfc(a,b,fc,fs);
% plot frequency response
a=anfc;
b=bnfc;
zz=exp(j*(0:fs/2-1)*2*pi/fs); % z values along unit circle
H=(b(1)*zz.^2+b(2)*zz+b(3))./(a(1)*zz.^2+a(2)*zz+a(3));
%H(H==0)=10^(-6); % avoid H=0 to avoid log of zero
H(1)=10^(-6);
figure(14);
subplot(2,1,1);
plot((1:fs/2)-1,abs(H));
axis([0 2*fc 0 max(abs(H))]);
subplot(2,1,2);
plot((1:fs/2),20*log10(abs(H/max(H))));
axis([0 fs/2 -40 0]);
title('|H(z)| 2-pole, 2-zero notch at fc');
xlabel('frequency (Hz)');
ylabel('amplitude');

% 2-pole 2-zero notch filter at 2*fc
% zeros at z1=e^+-j(theta),
% poles at p=|p|e^+-j(theta)
% where theta =  2*pi (2*fc)/fs
% for peak response at fc
% H(z) = (z-z1)(z-z1*)/(z-p)(z-p*)
%      = (z^2-(cos(theta)+1)/(z^2-2*|p|*cos(theta)+|p|^2)
theta6=2*pi*(2*fc)/fs;
p6=0.95; % sets bandwidth
% 2-sided 3dB BW(Hz) = (1-p)/(pi*sqrt(p))
gain6= (1-2*p6*cos(theta6)+p6^2)/(2-2*cos(theta6)); 
an2fc=[1 -2*p6*cos(theta6) p6*p6]; % denominator (poles)
bn2fc=gain6*[1 -2*cos(theta6) 1]; % numerator (zeros)
%plotBPFn2fc(an2fc,bn2fc,fc,fs);
%function plotBPFn2fc(a,b,fc,fs);
% plot frequency response
a=an2fc;
b=bn2fc;
zz=exp(j*(0:fs/2-1)*2*pi/fs); % z values along unit circle
H=(b(1)*zz.^2+b(2)*zz+b(3))./(a(1)*zz.^2+a(2)*zz+a(3));
%H(H==0)=10^(-6); % avoid H=0 to avoid log of zero
H(1)=10^(-6);
figure(15);
subplot(2,1,1);
plot((1:fs/2)-1,abs(H));
axis([0 4*fc 0 max(abs(H))]);
subplot(2,1,2);
plot((1:fs/2),20*log10(abs(H/max(H))));
axis([0 fs/2 -40 0]);
title('|H(z)| 2-pole, 2-zero notch at 2fc');
xlabel('frequency (Hz)');
ylabel('amplitude');

% ========================================
% channel
% ========================================
r_t=s_t; % no noise, filter or interference
% r_t=filter(bnfc,anfc,s_t); % filter out fc
% general channel
% r_t=g_t.*s_t + n_t + i_t
% g_t fading (time-varying gain)
% n_t noise, e.g. n_t=filter(b,a,randn(1,fs*T));
% i_t interference, e.g. a different s_t

% ========================================
% RECEIVER (demodulation)
% ========================================

% want to demodulate s_t to get mr_t which is
% the receiver's best estimate of m_t

% generate cos and sin carriers at receiver
phi=45*pi/180;  % phase offset
% phi=0;
% select ferr (difference between fc and
% carrier at receiver, ideally ferr=0)
% fixed value or ramp from ferrmin to ferrmax
% if want ramp
ferrmin=0;
ferrmax=8;
%ferr=ferrmin+(1:fs*T)*(ferrmax-ferrmin)/(fs*T);
% if want fixed value
ferr=0; % suggested values: 0,8,16
% frequency of cos and sin carriers
fcdel=fc+ferr;
% 2 pi fcdel t is written as below
twopi_fcdel_t=(1:fs*T)*2*pi.*fcdel/fs; 
cs_t=cos(twopi_fcdel_t+phi); % cos oscillator
sn_t=sin(twopi_fcdel_t+phi); % sin oscillator

%AM receiver
% rectifier, filter
%r_t(r_t<0)=0; % rectifier
%rf_t=filter(bfmax2,afmax2,r_t); % LPF cutoff fmax
%mr_t=rf_t; % estimate of m_t

% DSB receiver
%xr_t=r_t.*cs_t; % multiply with cs carrier
%xrf_t=filter(bfmax2,afmax2,xr_t); %LPF cutoff fmax
%mr_t=xrf_t; % estimate of m_t

% general i-q receiver
% obtains i,q,magnitude,phase
% an estimate of the modulation mr_t
% is obtained via the mapping for the 
% particular modulation type
xr_t=r_t.*cs_t; % i branch
% filter out double frequency term
xrf_t=filter(bfmax2,afmax2,xr_t); % LPF cutoff fmax
yr_t=r_t.*sn_t; % q branch
% filter out double frequency term
yrf_t=filter(bfmax2,afmax2,yr_t); % LPF cutoff fmax
% magnitude (envelope) of r_t
ar_t= xrf_t.^2+yrf_t.^2;
% phase of r_t
pr_t=atan2(yrf_t,xrf_t);
% for AM
%mr_t=ar_t; % and subtract DC component
% for DSB, SSB
mr_t=yrf_t;
% DSB frequency error signal 
% after LPF 5 Hz will be a DC level proportional
% to ferr
% e_t=xrf_t.*yrf_t;
% for FM
% mr_t = [diff(pr_t) 0]*fs/(2*pi); 
% why so ?
% instantaneous frequency 
% f_i(t)= (1/2*pi)*d(theta_i)/dt 
% theta_i(t) = 2 pi fc t + pr(t)
% f_i(t) = fc + (1/(2*pi))*d(pr_t)/dt
%        = fc + mr_t
% mr_t = (1/(2*pi))*d(pr_t)/dt
% approximate d(pr_t)=pr_t(i)-pr_t(i-1)
% approximate dt = t(i)-t(i-1) = ts=1/fs
% thus mr_t = [diff(pr_t) 0]*fs/(2*pi);
% add 0 to preserve length of mr_t

% for PM
% mr_t = pr_t;

% FM receiver (limiter, BPF, discriminator)
% limiter
% r_t(r_t<0)=-1;
% r_t(r_t>0)=1;
% BPF at fc
% r_t=filter(bfc,afc,r_t); % BPF at fc
% choose one detection method
% one of discriminator, delay, slope 
    % discriminator (differentiator)
       % r_t=[diff(r_t) 0]*fs;
    % delay 0.5 of a carrier wavelength
       % delay = zeros(1,(fs/(2*fc))-2); 
       % method 1:
       % dfilter=[1 delay -1];
       % r_t=filter(dfilter,1,r_t);
        % method 2: 
        % delayed r_t (same length as r_t)
        %r_tdelay = [delay  r_t(1:(length(r_t)-length(delay)))];
        % r_t = r_t - r_tdelay;
    % slope detection on lower slope of BPF at 2fc
        % r_t=filter(b2fc,a2fc,r_t);
% now have AM signal
%r_t(r_t<0)=0; % rectify
%rf_t=filter(bfmax2,afmax2,r_t); %LPF cutoff fmax
%mr_t=rf_t; % estimate of m_t + DC term

% ====================================
% FFT
% ====================================

% take fft of s_t and mr_t
% first divide s_t into blocks of length blockfft
ys=reshape(s_t,blockfft,fs*T/blockfft);
yr=reshape(xrf_t,blockfft,fs*T/blockfft);
yf=reshape(mr_t,blockfft,fs*T/blockfft);
ym=reshape(m_t,blockfft,fs*T/blockfft);
%yd=reshape(cdif_t,blockfft,fs*T/blockfft);
yc=reshape(c_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(ys)/blockfft; % fft operation applied to each column
zm=fft(yf)/blockfft; % 

% ========================================
% generate plots (time and frequency domain)
% ========================================

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

figure(2); % 3-D view of FFTs of s_t
% set up xx and yy axes for 3-d mesh plot
xx=(0:blockfft-1)*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); % view(0,90) for top view
title('3-D view of FFTs of s(t)');
xlabel('frequency (Hz)');
zlabel('amplitude');
ylabel('block #');

figure(11); % 2-D view of FFTs of mr_t
zma=zm(:);maxzm=max(abs(zma));% find max of all ffts
%plot((1:blockfft)*fs*T/blockfft,abs(z(:,1)));
plot((0:blockfft-1)*fs/blockfft,abs(zm));
axis([0 2*fmax 0 maxzm]);
title('2-D view of FFTs of mr(t)');
xlabel('frequency (Hz)');
ylabel('amplitude');

figure(12); % 3-D view of FFTs of mr_t
% set up xx and yy axes for 3-d mesh plot
xx=(0:blockfft-1)*fs/blockfft;
yy=(1:fs*T/blockfft);
zmt=zm';
mesh(xx,yy,abs(zmt)); % makes 3-d plot
axis([0 2*fmax 0 fs*T/blockfft 0 maxzm]);
view(25,30); %view(0,90) for top view
title('3-D view of FFTs of mr(t)');
xlabel('frequency (Hz)');
zlabel('amplitude');
ylabel('block #');


figure(3); %time domain plot of first block
plot((1:blockfft)/fs,ys(:,1),'b-');hold on
% if no receiver, comment out yr, yf
plot((1:blockfft)/fs,yr(:,1),'g-');
plot((1:blockfft)/fs,yf(:,1),'m-');
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
title('middle block (time domain)');
xlabel('time (sec)');
ylabel('amplitude');
subplot(2,2,1);
plot((1:blockfft)/fs,ys(:,nblocks/2),'b-');
title('s(t)');
xlabel('time (sec)');
ylabel('amplitude');
subplot(2,2,2);
plot((1:blockfft)/fs,yr(:,nblocks/2),'g-');
title('xr(t)');
xlabel('time (sec)');
ylabel('amplitude');
subplot(2,2,3);
plot((1:blockfft)/fs,yf(:,nblocks/2),'m-');
title('mr(t)');
xlabel('time (sec)');
ylabel('amplitude');
subplot(2,2,4);
plot((1:blockfft)/fs,ym(:,nblocks/2),'r-');
title('m(t)');
xlabel('time (sec)');
ylabel('amplitude');


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

% output sound file of s_t and mr_t
s_tmax=max(abs(s_t));
%soundsc(s_t/s_tmax,fs);
%wavwrite(s_t/s_tmax,fs,'amfmdem3.wav'); % for win
mr_tmax=max(abs(mr_t));
%wavwrite(mr_t/mr_tmax,fs,'amfm3mr.wav'); % for win

%auwrite(s_t/s_tmax,fs,'amfm5.au'); % for unix

%testint(twopi_fmod_t,fs,fmod);
% tests numerical integration and differentiation
%function testint(twopi_fmod_t,fs,fmod);
% numerical integration and differentiation
% test integral and derivative of cos(twopi_f_t)
% and shows they are equal to sin(twopi_f_t)
% for sampled signals with sampling time ts=1/fs
% the infinistesimal dt becomes ts=1/fs
% int f_t dt = cumsum(f_t) * ts = cumsum(f_t)/fs
% deriv f_t/dt = diff(f_t) / ts = diff(f_t)*fs
% only works for fixed scalar fmod
% fcn = sin(twopi_fmod_t);
% int = 2*pi*fmod*cumsum(cos(twopi_fmod_t)/fs);
% deriv = [diff(-cos(twopi_fmod_t)) 0]*fs/(2*pi*fmod);
% figure(13);
% plot(fcn(1:1024),'b-'); hold on;
% plot(int(1:1024),'r-'); 
% plot(deriv(1:1024),'g-');
% hold off;

keyboard; % preserve interactive mode
% so can look at value of variables from
% command line
bcd=4; % dummy command;





--=====================_908811974==_--