%  This script demonstrates the modulation properties of the real and
%  complex sinusoid.  A simple sinusoid is used as a test signal to allow
%  for easy observation of the frequency shifts involved in the modulation
%  process.  Alternatively, certain lines can be uncommented to read in
%  a speech wavefile and modulate that (instead of the simple sinusoid) and
%  observe the spectral shifts.  The sampling rate should be set
%  sufficiently high to account for the shift of the modulation frequency
%  plus the frequency range of the original signal.  The PSD is estimated
%  with the PWELCH function to observe the spectrum (the average spectrum
%  better for plot when signal other than the sinusoid are used) and
%  therefore a hopping window length is selected.  If the signal is
%  deterministic (i.e. a pulse or short segment waveform), the window
%  length should be set equal to the signal duration.  In this case PWELCH
%  will just compute the FFT magnitude (no averaging over multple windows
%  will occur).  The 'twosided' option is selected for PWLECH and the 
%  FFTSHIFT is used to see both the positive and negative frequencies on
%  an axis symmetric about zero.  This is important to observe the shifting
%  of both the positive and negative frequency components as well as
%  handling the complex valued waveforms, whose spectra are not symmetric
%  about 0 (DC).
%  This script generates plots of the spectra magnitudes for test signal,
%  modulating (carrier) frequency, and the modulated signal for the real
% (symetric spectra) and complex-valued (one-sided spectra) signals.
%
%

%  Set parameters
fs = 2*48e3;  %  Sampling frequency in Hz (signal read in will be resampled)
ft = 4e3;   %  Test frequency a Sinudoid in Hz (applies if sin signal code is uncommented)
fc = 15000;  %  Carrier Frequency a Sinusoid in Hz
tdur = 1;  %  Duration of signal in seconds
wl = round(100*fs/ft);  %  Window size for PSD contains 20 cycles of test frequency 

%  Generate real signal example



%  Test Waveform   %  UNCOMMENT NEXT 4 LINES TO READ IN SIGNAL
% [st, fso] = wavread('ImpNtrans.wav');    
% st = st'; 
% st = resample(st,fs,fso);  %  Resample signal to higer sampling rate
% tax = (0:length(st)-1)/fs;   % Time axis

% Simple sinusoid test signal  %   UNCOMMENT NEXT 2 LINES TO USE SINUSOID 
tax = (0:round(tdur*fs)-1)/fs;  % Time axis
st = sin(2*pi*ft*tax);   

%  Carrier Waveform
sc = cos(2*pi*fc*tax);

%  Modulation
modsig = sc.*st;


%  Compute Spectra (two-sided)
[psc, f] = pwelch(sc,hamming(wl),round(wl/2), 2*wl, fs, 'twosided');
[pst, f] = pwelch(st,hamming(wl),round(wl/2), 2*wl, fs, 'twosided');
[pmodsig, f] = pwelch(modsig,hamming(wl),round(wl/2), 2*wl, fs, 'twosided');
%  Create frequency axis symetric about 0 
fsym = fs*[-round(length(f)/2):round(length(f)/2)-1]/length(f);
%  Plot Spectra using FFTSHIFT to make symetric about 0
figure(1)
subplot(3,1,1); plot(fsym/1000,10*log10(fftshift(psc)*fs/wl));
xlabel('kHz'); ylabel('dB'); title('Carrier Frequency')
set(gca,'Ylim', [-90 0],'Xlim', [-fs/2 fs/2]/1000);
subplot(3,1,2); plot(fsym/1000,10*log10(fftshift(pst)*fs/wl)); 
xlabel('kHz'); ylabel('dB'); title('Test Signal')
set(gca,'Ylim', [-90 0],'Xlim', [-fs/2 fs/2]/1000);
subplot(3,1,3); plot(fsym/1000,10*log10(fftshift(pmodsig)*fs/wl))
xlabel('kHz'); ylabel('dB'); title('Modulated Signal')
set(gca,'Ylim', [-90 0],'Xlim', [-fs/2 fs/2]/1000);


%  Generate Complex Signal example 

%  Complex Carrier Waveform
scx = exp(j*2*pi*fc*tax);
%  make test signal complex so that spectrum is single sided
stx = hilbert(st);
%  Perform complex modulation 
modsigx = scx.*stx;

%  Compute Spectrum (two-sided)
[pscx, f] = pwelch(scx,hamming(wl),round(wl/2), 2*wl, fs, 'twosided');
[pstx, f] = pwelch(stx,hamming(wl),round(wl/2), 2*wl, fs, 'twosided');
[pmodsigx, f] = pwelch(modsigx,hamming(wl),round(wl/2), 2*wl, fs, 'twosided');
%  Create frequency axis symetric about 0 
fsym = fs*[-round(length(f)/2):round(length(f)/2)-1]/length(f);
%  Plot Spectra using FFTSHIFT to make symetric about 0
figure(2)
subplot(3,1,1); plot(fsym/1000,10*log10(fftshift(pscx)*fs/wl));
xlabel('kHz'); ylabel('dB'); title('Carrier Frequency')
set(gca,'Ylim', [-90 0],'Xlim', [-fs/2 fs/2]/1000);
subplot(3,1,2); plot(fsym/1000,10*log10(fftshift(pstx)*fs/wl)); 
xlabel('kHz'); ylabel('dB'); title('Test Signal')
set(gca,'Ylim', [-90 0],'Xlim', [-fs/2 fs/2]/1000);
subplot(3,1,3); plot(fsym/1000,10*log10(fftshift(pmodsigx)*fs/wl))
xlabel('kHz'); ylabel('dB'); title('Modulated Signal')
set(gca,'Ylim', [-90 0],'Xlim', [-fs/2 fs/2]/1000);
