Updating documentation and adding examples.

doc/make.jl

@@ -1,5 +1,17 @@
 push!(LOAD_PATH, "../src/")
 using Documenter, DigitalComm
 
-makedocs(sitename="DigitalComm")
+makedocs(sitename="DigitalComm.jl", 
+		 format = Documenter.HTML(prettyurls = false),
+		 pages    = Any[
+						"Introduction to DigitalComm"   => "index.md",
+						"Function list"          => "base.md",
+						"Examples"                => Any[ 
+														 "Examples/example_AWGN.md",
+														 "Examples/example_BER.md",
+														 "Examples/example_PSD.md",
+														 ],
+						],
+		 );
+
 #makedocs(sitename="My Documentation", format = Documenter.HTML(prettyurls = false))

doc/src/Examples/example_AWGN.md

@@ -0,0 +1,44 @@
+# Transmission of xQAM with additive white Gaussian noise 
+
+To simulate a transmission of QPSK // 16QAM // 64QAM // 256QAM other a
+white additive Gaussian noise and display the received constellation,
+the following code can be used 
+
+    # ---------------------------------------------------- 
+    # --- Transmitter
+    # ---------------------------------------------------- 
+    using DigitalComm 
+    using Plots 
+    # --- Parameters 
+    snr		  = 20;
+    mcs		  = 16;
+    nbBits	  = 1024* Int(log2(mcs));
+    # --- Binary sequence generation 
+    bitSeq	  = genBitSequence(nbBits);
+    # --- QPSK mapping
+    qamSeq	  = bitMappingQAM(mcs,bitSeq);
+    # ---------------------------------------------------- 
+    # --- Channel  
+    # ---------------------------------------------------- 
+    #  --- AWGN
+    # Theoretical power is 1 (normalized constellation)
+    qamNoise,  = addNoise(qamSeq,snr,1);
+    # ----------------------------------------------------
+    # --- Rx Stage: SRRC
+    # ----------------------------------------------------
+    # --- Binary demapper
+    bitDec	= bitDemappingQAM(mcs,qamNoise);
+    # --- BER measure
+    ber	  = sum(xor.(bitDec,bitSeq)) /length(bitSeq);
+    # --- Display constellation 
+    plt	  = scatter(real(qamNoise),imag(qamNoise),label="Noisy");
+	scatter!(plt,real(qamSeq),imag(qamSeq),label="Ideal");
+	xlabel!("Real part");
+	ylabel!("Imag part");
+    display(plt);
+
+It plots the received constellation impaired by noise (here a 20dB SNR is used) 
+![Constellation](./../img/constellation.png)
+
+
+

doc/src/Examples/example_BER.md

@@ -0,0 +1,121 @@
+# Compute the theoretical BER for AWGN channel and various constellation size 
+
+Based on the previous skeleton, we can now compute an iterative testbench to compute the Bit Error Rate for various constellation size, and compare the simulation with the theory.
+As  a gentle reminder, the theoretical bit error rate can be approximated as 
+
+``\mathrm{BER} = \frac{ 4 \left( 1 - \frac{1}{\sqrt{M}} \right) }{ \log_2(M)} \times Q( \sqrt{ \frac{6  \log_2(M)}{2(M-1)} \frac{Eb}{N_0}}``
+
+
+First of all let's call the modules 
+
+
+	using DigitalComm 
+	using PGFPlotsX
+
+
+We define first the main monte-carlo function that compute an elementary Tx-Rx link, and returns the number of error and number of bit computed (to be accumulated) 
+
+
+	function monteCarlo(snr,mcs,nbSymb) 
+	  # Number of bits  
+	  nbBits	  = nbSymb * Int(log2(mcs));
+	  # --- Binary sequence generation 
+	  bitSeq	  = genBitSequence(nbBits);
+	  # --- QPSK mapping
+	  qamSeq	  = bitMappingQAM(mcs,bitSeq);
+	  # ---------------------------------------------------- 
+	  # --- Channel  
+	  # ---------------------------------------------------- 
+	  #  --- AWGN
+	  # Theoretical power is 1 (normalized constellation)
+	  qamNoise,  = addNoise(qamSeq,snr,1);
+	  # ----------------------------------------------------
+	  # --- Rx Stage: SRRC
+	  # ----------------------------------------------------
+	  # --- Binary demapper
+	  bitDec	= bitDemappingQAM(mcs,qamNoise);
+	  # --- Error counter 
+	  nbE		= sum(xor.(bitDec,bitSeq));
+	  # --- Return Error and bits 
+	  return (nbE,nbBits);
+	end
+
+
+A function to plot the BER versus the SNR, for different mcs and compare to theory 
+
+
+	function doPlot(snrVect,ber,qamVect)
+		a = 0;
+		@pgf a = Axis({
+					   ymode	  = "log",
+					   height      ="3in",
+					   width       ="4in",
+					   grid,
+					   xlabel      = "SNR [dB]",
+					   ylabel      = "Bit Error Rate ",
+					   ymax 	   = 1,
+					   ymin 	   = 10.0^(-5),
+					   title       = "AWGN BER for QAM",
+					   legend_style="{at={(0,0)},anchor=south west,legend cell align=left,align=left,draw=white!15!black}"
+					   },
+					  Plot({color="red",mark="square*"},Table([snrVect,ber[1,:]])),
+					  LegendEntry("QPSK"),
+					  Plot({color="green",mark="*"},Table([snrVect,ber[2,:]])),
+					  LegendEntry("16-QAM"),
+
+					  Plot({color="purple",mark="triangle*"},Table([snrVect,ber[3,:]])),
+					  LegendEntry("64-QAM"),
+					  Plot({color="blue",mark="diamond*"},Table([snrVect,ber[4,:]])),
+					  LegendEntry("256-QAM"),
+					  );
+		# ---  Adding theoretical curve
+		snrLin  = (10.0).^(snrVect/10)
+		for qamScheme = qamVect
+			ebNo 	= snrLin / log2(qamScheme);
+			# This approximation is only valid for high SNR (one symbol error is converted to one bit error with Gray coding).
+			berTheo	  = 4 * ( 1 - 1 / sqrt(qamScheme)) / log2(qamScheme) * qFunc.(sqrt.( 2*ebNo * 3 * log2(qamScheme) / (2*(qamScheme-1)  )));
+			@pgf push!(a,Plot({color="black"},Table([snrVect,berTheo])));
+		end
+		display(a);
+	end 
+
+
+Then, the main routine to compute the BER for a given number of iterations and a range of SNR 
+
+
+	function main() 
+	# --- Parameters 
+	nbIt		    = 10000;			  # Number of iterations  
+	nbSymb	        = 1024;			  # Number of symbols per iterations 
+	mcs		        = [4,16,64,256];  # Constellation size 
+	snrRange	    = (-1:26);		  # SNR, expressed in dB 
+	# --- Init performance metrics
+	nbSNR			= length(snrRange);
+	ber			    = zeros(Float64,length(mcs),nbSNR);
+	for iMcs = 1 : 1  : length(mcs)
+		for iSNR = 1 : 1 : nbSNR
+			# --- Create BER counters 
+			nbE	  = 0;
+			nbB	  = 0;
+			for iN = 1 : 1 : nbIt
+				# --- Elementary MC call 
+				# Corresponds to a given SNR and a given iteration 
+				# As we are ergodic in AWGN, it is only nbSymb*nbIt that matters for BER computation
+				(a,b) = monteCarlo(snrRange[iSNR],mcs[iMcs],nbSymb); 
+				# --- Update counters 
+				nbE	+= a;	  # Increment errors 
+				nbB	+= b;	  # Increment bit counters 
+			end 
+			ber[iMcs,iSNR] = nbE / nbB; 
+		end
+	end 
+	# --- Plotting routine
+	doPlot(snrRange,ber,mcs);
+	end
+	  
+
+The output plot is the following, showing adequacy between theory and practise for high SNR (the theoretical curve is under the assumption that one symbol error leads to one erroneous bit (gray coding) which is true only with intermediate noise levels).
+![BER](./../img/BER_AWGN.png)
+
+
+

doc/src/Examples/example_PSD.md

@@ -0,0 +1,184 @@
+# Plotting PSD of several multicarrier waveforms 
+
+The purpose is to compre the Power spectral density of several multicarrier waveform. The following module can be used: 
+
+	module example_PSD_waveform 
+	# ---------------------------------------------------- 
+	# --- Modules  
+	# ---------------------------------------------------- 
+	using DigitalComm 
+	# --- External Modules
+	using Plots 
+	gr();
+	using Printf
+	using FFTW
+	# ---------------------------------------------------- 
+	# --- Core functions  
+	# ---------------------------------------------------- 
+	""" psdWaveform.m 
+	---  
+	Compute the power spectral density (i.e the spectrum here) of the signal parametrized by the waveform structure waveform, for a number of symbol nbSymb.
+	The frequency allocation is the one inherited from the waveform structure (i.e waveform.allocatedSubcarriers).
+	# --- Syntax 
+	( freq,psd )	= psdWaveform(waveform,nbSymb,allocatedSubcarriers);
+	# ---  Input parameters 
+	- waveform  : Structure associated to transmitted waveform 
+	- nbSymb	  : Number of symbol to be transmitted [Int]
+	- nbIt	  : Monte carlo parameter for PSD evaluation (should be > 1)
+	# --- Output parameters 
+	- freq	  : Vector of frequency evaluation (between -0.5 and 0.5). [Array{Float64,L}]
+	- psd		  : Spectrum evaluated on freq [Array{Complex{Float64}},L]
+	# --- Input parameters 
+	- 
+	# --- Output parameters 
+	- 
+	# --- 
+	# v 1.0 - Robin Gerzaguet.
+	"""
+	function  psdWaveform(waveform,nbSymb,nbIt)
+		# ----------------------------------------------------
+		# --- PSD calculation
+		# ---------------------------------------------------- 
+		# --- Getting frequency allocation 
+		allocatedSubcarriers  = waveform.allocatedSubcarriers;
+		# --- Getting number of bits 
+		# First, frequency size 
+		nbSubcarriers	      = length(allocatedSubcarriers);
+		# Force a fiven mcs 
+		mcs				      = 4;	  # QPSK.
+		# Deduce number of required bits 
+		nbBits			      = nbSymb * nbSubcarriers * Int(log2(mcs));
+		# --- Init psd evaluator 
+		psd = 0;
+		# --- Iterative PSD calculation
+		for iN = 1 : 1 : nbIt
+			# --- Binary sequence
+			bitSeq	      = genBitSequence(nbBits);
+			# Mapping
+			qamSeq		  = bitMappingQAM(mcs,bitSeq);
+			# --- T/F matrix
+			qamMat		  = reshape(qamSeq,nbSubcarriers,nbSymb);
+			# --- Signal
+			sigPSD		  = genSig(qamMat,waveform);
+			# --- Mean PSD:
+			psd		  	  = psd .+ 1/nbIt*1/length(sigPSD)*abs.(fftshift(fft(sigPSD))).^2;
+		end
+		# --- Calculating sampling frequency
+		# Returns Nyquist frequency 
+		fe		= 1;
+		Basefe	= (0:(length(psd) .-1))./length(psd)*fe .-fe/2;
+		return (Basefe,psd);
+	end
+	# ---------------------------------------------------- 
+	# --- Main routine  
+	# ---------------------------------------------------- 
+	function main()
+		# ----------------------------------------------------
+		# --- Overall parameters
+		# ----------------------------------------------------
+		# --- Overall PHY parameters
+		nbIt			= 50;			  # --- Iteration number
+		nbSymb 			= 14;			  # --- Number of symbols (one frame)
+		nFFT 			= 1024;			  # --- Base FFT size
+		samplingFreq	= 15.36;		  # --- Frequency value (MHz)
+		# --- Frequency allocation
+		#allocatedSubcarriers= getLTEAlloc(nFFT);
+		#allocatedSubcarriers = (1:12*4);
+		# 4 RB alloc. 1 RB space. 4 RB allocated
+		allocatedSubcarriers = [1:12*4; 12*5 .+ (1:12*4)];
+		# ----------------------------------------------------
+		# --- Waveform contender
+		# ----------------------------------------------------
+		# --- Init OFDM structure
+		ofdm  = initOFDM(
+						 nFFT,						        # --- nFFT                 : FFT size
+						 72,						        # --- nCP                  : CP size
+						 allocatedSubcarriers		        # --- allocatedSubcarriers : Subcarrier allocation
+						 );
+		# --- Init SCFDMA structure
+		scfdma  = initSCFDMA(
+							 nFFT,						        # --- nFFT                 : FFT size
+							 72,						        # --- nCP                  : CP size
+							 allocatedSubcarriers,		        # --- allocatedSubcarriers : Subcarrier allocation
+							 12;								# --- sizeDFT			   : DFT preprocessing size
+							 );
+		# --- Init UF-OFDM structure
+		ufofdm  = initUFOFDM(
+							 nFFT,					        # --- nFFT                 : FFT size
+							 73,					        # --- L                    : Filter length (same size +1 due to conv)
+							 allocatedSubcarriers,	        # --- allocatedSubcarriers : Subcarrier allocation
+							 applyPD=1,				        # --- applyPD              : Do predistortion at Tx stage
+							 attenuation=40,		        # --- attenuation          : Filter attenuation in dB
+							 );
+		# --- Init BF-OFDM structure
+		bfofdm	= initBFOFDM(
+							 32,				            # --- nFBMC                : PPN size (max number of carriers)
+							 64,				            # --- nOFDM                : Precoder size (OFDM sizer)
+							 3,				            	# --- K                    : Overlapping factor
+							 9,					            # --- GI                   : CP size of precoder
+							 0.5,				            # --- δ                    : compression factor
+							 allocatedSubcarriers,          # --- allocatedSubcarriers : Subcarrier allocation
+							 "gaussian",	            # --- filterName           : Pulse shape name
+							 BT=0.36,				        # --- BT                   : Potential BT value for Gaussian
+							 filterStopBand = 110,			# --- filterStopBand       : DC stopband value
+							 fS=[],				            # --- fS                   : Potential frequency coefficient for FS filter
+							 nFFT= 1024,		            # --- nFFT                 : associated FFT value in Rx
+							 nCP= 72,			            # --- nCP                  : extended CP size
+							 );
+		# --- Init WOLA-OFDM structure
+		wola  = initWOLA(
+						 nFFT,						        # --- nFFT                 : FFT size
+						 72,						        # --- nCP                  : CP size
+						 allocatedSubcarriers,		        # --- allocatedSubcarriers : Subcarrier allocation
+						 "triangle",						# --- Window type @Tx side
+						 20,								# --- Window size @Tx side
+						 "triangle",						# --- Window type @Rx side
+						 20,								# --- Window size @Rx side
+						 );
+		fbmc  = initFBMC(
+						 nFFT,						        # --- nFFT                 : FFT size
+						 4,									# --- K					   : Overlapping factor
+						 allocatedSubcarriers		        # --- allocatedSubcarriers : Subcarrier allocation
+						 );
+		# ----------------------------------------------------
+		# --- Merging structures
+		# ----------------------------------------------------
+		# Create  a dictionnary to rule them all 
+		waveforms 	= initWaveforms(ofdm,
+									scfdma,
+									ufofdm,
+									bfofdm,
+									wola,
+									fbmc,
+									);
+		# ---------------------------------------------------- 
+		# --- PSD main calculation  
+		# ---------------------------------------------------- 
+		# --- Init plot container 
+		plt	  = plot(reuse=false);
+		decim	= 1;	# decimation for light plots
+		# --- Iterative PSD generation
+		for (name,struc) in waveforms 
+			# --- Calculate PSD for the configuration 
+			(fe,psd)  = psdWaveform(struc,nbSymb,nbIt); 
+			# Plot the result 
+			plot!(plt,fe[1:decim:end].*samplingFreq,10 .* log10.(psd[1:decim:end]/maximum(psd)),label=name,legend=:topleft);
+		end
+		# --- Update plot and adding labels 
+		# Purpose is to zoom out on allocated region.
+		scsN 		= (1/1024)*samplingFreq; 		# Subscarrier spacing (normalized)
+		rbV 		= (12*12); 		# See several RB for psd fall-off
+		ylims!(-120,5);
+		xlims!(-rbV*scsN,maximum(allocatedSubcarriers)*scsN+2*12*scsN);
+		xlabel!("Frequency [MHz]");
+		ylabel!("Spectrum");
+		display(plt)
+	end
+	end
+
+By running `example_PSD_waveform.main();` a comparison plot between the different PSD can be obtained 
+
+![PSD](./../img/psd.png)
+
+
+