Adding CDMA to signal generation

examples/example_BER_CDMA.jl

@@ -0,0 +1,151 @@
+module example_BER_cdma
+
+
+
+# ---------------------------------------------------- 
+# --- Modules  
+# ---------------------------------------------------- 
+using DigitalComm 
+# --- External Modules
+using Plots
+using Printf
+using FFTW
+using Statistics
+
+
+
+# ---------------------------------------------------- 
+# --- Core functions  
+# ---------------------------------------------------- 
+function  getBER(waveform,mcs,snrVect,nbSymb,nbIt)
+    if waveform isa DigitalComm.StrucCDMA 
+        # --- Get code allocation 
+        nbSubcarriers       = length(waveform.userMask)
+    else 
+        # --- Getting frequency allocation 
+        allocatedSubcarriers  = waveform.allocatedSubcarriers;
+        nbSubcarriers	      = length(allocatedSubcarriers);
+    end
+	# Deduce number of required bits 
+	nbBits			      = nbSymb * nbSubcarriers * Int(log2(mcs));
+	# --- Init BER vector
+	nbSNR			= length(snrVect);
+	ber				= zeros(Float64,nbSNR);
+	# --- Iterative PSD calculation
+	for k = 1 : 1 : nbSNR
+		# --- Update counters
+		nbC		= 0;
+		nbE		= 0;
+		powSig	= 0;
+		for iN = 1 : 1 : nbIt
+			# --- Binary sequence
+			bitSeq	      = genBitSequence(nbBits);
+			# Mapping
+			qamSeq		  = bitMappingQAM(mcs,bitSeq);
+			# --- T/F matrix
+			qamMat		  = reshape(qamSeq,nbSubcarriers,nbSymb);
+			# --- Signal
+			sigId		  = genSig(qamMat,waveform);
+			# ---------------------------------------------------- 
+			# --- Channel  
+			# ---------------------------------------------------- 
+			#  --- AWGN
+			if iN == 1 
+				# We compute the power based on generated signal 
+				# We troncate the beginning and end of signal to avoid 
+				# estimation biais induced by (potential) filter tails
+				powSig = mean(abs2.( @views sigId[1+ end÷4 : end - end÷4]));
+			end
+			sigNoise,  = addNoise(sigId,snrVect[k],powSig);
+			# ----------------------------------------------------
+			# --- Rx Stage
+			# ----------------------------------------------------
+			# --- Waveform demodulator 
+			qamDec	  = decodeSig(sigNoise,waveform);
+			# --- Binary demapper
+			bitDec	  = bitDemappingQAM(mcs,qamDec[:]);
+			# --- BER measure
+			nbE	 += sum(xor.(bitDec,bitSeq));
+			nbC	 += length(bitSeq);
+		end
+		# --- BER measure
+		# Adding 1e-10 to avoid log plot of zero errors
+		ber[k]		= 1e-10 .+ nbE / nbC;
+	end
+	return ber
+end
+
+# ---------------------------------------------------- 
+# --- Main routine  
+# ---------------------------------------------------- 
+function main()
+	# ----------------------------------------------------
+	# --- Overall parameters
+	# ----------------------------------------------------
+	# --- Overall PHY parameters
+	nbIt			= 20;			  # --- Iteration number
+	nbSymb 			= 1400;			  # --- Number of symbols (one frame)
+	nFFT 			= 1024;			  # --- Base FFT size
+	samplingFreq	= 3.84e6;		  # --- Frequency value (MHz)
+	mcs				= 16;			  # --- 16-QAM 
+	snrVect			= (-10:30);
+	# --- Frequency allocation
+	allocatedSubcarriers= getLTEAlloc(nFFT);
+
+	# ----------------------------------------------------
+	# --- Waveform contender
+	# ----------------------------------------------------
+	# --- Init OFDM structure
+    cdma4 = initCDMA(
+                    4,
+                    :ovsf
+                   )
+    cdma16 = initCDMA(
+                    16,
+                    :ovsf
+                   )
+    cdma32 = initCDMA(
+                    32,
+                    :ovsf
+                   )
+    cdma64 = initCDMA(
+                    64,
+                    :ovsf
+                   )
+	# ----------------------------------------------------
+	# --- Merging structures
+	# ----------------------------------------------------
+	# Create  a dictionnary to rule them all 
+	waveforms 	= initWaveforms(
+                                cdma4,
+                                cdma16,
+                                cdma32,
+                                cdma64,
+								);
+
+
+	# ---------------------------------------------------- 
+	# --- BER main calculation  
+	# ---------------------------------------------------- 
+	# --- Init plot container 
+	plt	  = plot(reuse=false,yscale=:log10,legend=:bottomleft);
+	# --- Iterative PSD generation
+	for (name,struc) in waveforms 
+		# --- Calculate PSD for the configuration 
+		ber	= getBER(struc,mcs,snrVect,nbSymb,nbIt);
+		# --- Plot stuff
+		plot!(plt,snrVect,ber,label=struc.nbUsers);
+	end
+	# --- Adding metata do plot curve
+	xlabel!("SNR [dB]");
+	ylabel!("Bit Error Rate");
+	ylims!(1e-6,1);
+	display(plt);
+end
+
+
+
+end
+
+
+

examples/example_cdma.jl

@@ -0,0 +1,42 @@
+using DigitalComm 
+
+
+# ----------------------------------------------------
+# --- Parameters 
+# ---------------------------------------------------- 
+nbBits  = 32768 
+mcs     = 4 
+nbUsers = 16
+nbActiveUsers = 16
+
+# ----------------------------------------------------
+# --- CDMA Generation
+# ---------------------------------------------------- 
+# --- Binary sequence
+bitSeq	      = genBitSequence(nbBits)
+# Mapping
+qamSeq		  = bitMappingQAM(mcs,bitSeq);
+nbSymb        = length(qamSeq) ÷ nbActiveUsers
+# --- T/F matrix
+qamMat		  = reshape(qamSeq,nbActiveUsers,nbSymb);
+# --- Signal
+sigId		  = cdmaSigGen(qamMat,nbUsers,:ovsf)
+
+
+# ----------------------------------------------------
+# --- Ideal decoding
+# ---------------------------------------------------- 
+# --- Decoding 
+qamRx = cdmaSigDecode(sigId,nbUsers,:ovsf,1:nbActiveUsers)
+sir = getSIR(qamRx,qamMat)
+println("SIR between Tx and Rx sequence (no noise) is $sir dB")
+
+
+# ----------------------------------------------------
+# --- Decoding with noise
+# ---------------------------------------------------- 
+snr = 30 
+sigRx,_ = addNoise(sigId,snr)
+qamRx = cdmaSigDecode(sigRx,nbUsers,:ovsf,1:nbActiveUsers)
+sir = getSIR(qamRx,qamMat)
+println("SIR between Tx and Rx sequence ($snr dB additive  noise) is $sir dB")

src/DigitalComm.jl

@@ -64,20 +64,15 @@ export sqrtRaisedCosine
 # ---------------------------------------------------- 
 # --- Q function definition 
 import SpecialFunctions.erfc
-"""
----  
-Returns the Q function used for Bit error rate computation in digital system 
- Q(x)= 1/2 erfc(x/sqrt(2))
-erfc is the Complexmplementary error function, i.e. the accurate version of 1-erf(x) for large x
-erfc is inherited from DSP
-# --- Syntax 
-y = qfunc(x)
-# --- Input parameters 
+"""Returns the Q function used for Bit error rate computation in digital system \\
+ Q(x)= 1/2 erfc(x/sqrt(2)) \\
+erfc is the Complexmplementary error function, i.e. the accurate version of 1-erf(x) for large x \\
+erfc is inherited from DSP\\
+y = qfunc(x) \\
+Input parameters 
 - x: Input [Float64]
-# --- Output parameters 
+Output parameters 
 - y: Q(x)[Float64] 
-# --- 
-# v 1.0 - Robin Gerzaguet.
 """
 function qFunc(x) 
 	return 1/2 * erfc.( x / sqrt(2));
@@ -87,21 +82,18 @@ export qFunc;
 
 
 """
---- 
-Returns the Signal to interference ratio expressed in dB (or in linear) between a obersvation signal d(n) and a reference signal u(n)
-The ratio is expressed as 10*log10( E[ || d(n) - u(n) || / E[||u(n)||^2]  )
- with E the expectation wrt to time
- The 2 vectors d and u should have the same length L
-# --- Syntax 
-  sir = getSIR( d, u , type="dB")
-# ---  Input parameter 
-- d	: Observation signal [Array{Any}]
-- u	: Reference signal [Array{Any}]	
-- type: Output unit [String]: "dB" or "Linear" (default, "dB")
-# --- Output parameters
+Returns the Signal to interference ratio expressed in dB (or in linear) between a obersvation signal d(n) and a reference signal u(n) \\
+The ratio is expressed as 10*log10( E[ || d(n) - u(n) || / E[||u(n)||^2]  ) \\
+ with E the expectation wrt to time \\
+The 2 vectors d and u should have the same length L \\
+
+sir = getSIR( d, u , type="dB")\\
+Input parameters 
+- d	: Observation signal [`Array`]
+- u	: Reference signal [`Array`]	
+- type: Output unit [`String`]: "dB" or "Linear" (default, "dB")
+Output parameters
 - sir	: Signal to interference ratio in unit `type`
-# --- 
-# v 1.0 - Robin Gerzaguet.
 """
 function getSIR(d,u,type="dB")
 	# --- Setting type
@@ -233,6 +225,14 @@ include("./Waveforms/FBMC/fbmcSigDecode.jl");
 export oqamDemapping 
 export fbmcSigDecode
 
+# CDMA Generation 
+include("./Waveforms/CDMA/cdmaSigGen.jl")
+export initCDMA
+export ovsf 
+export cdmaSigGen
+export cdmaSigGen!
+include("./Waveforms/CDMA/cdmaSigDecode.jl")
+export cdmaSigDecode
 
 
 # --- Global waveform alias

src/Waveforms/CDMA/cdmaSigDecode.jl

@@ -0,0 +1,38 @@
+
+
+function cdmaSigDecode(sig::AbstractVector{T},nbUsers,code::Symbol,userMask=nothing) where T
+    # --- Define constant 
+    nbChips = length(sig)
+    nbSymb = nbChips  ÷ nbUsers
+    # --- Define code 
+    if code == :ovsf
+        c = ovsf(nbUsers)
+    else 
+        @error "Code $code is not supported"
+    end
+    if isnothing(userMask)
+        userMask = 1:nbUsers
+    end
+    nbActiveUsers = length(userMask)
+    qamRx = zeros(T,nbActiveUsers,nbSymb)
+    tmp = zeros(T,nbSymb) # To be optimized as a view of qamRx
+    for n ∈ eachindex(userMask)
+        despreading!(tmp,sig,nbUsers,c[:,userMask[n]])
+        qamRx[n,:] .= tmp 
+    end
+    return qamRx 
+end
+
+function despreading!(tmp,sig,nbUsers,code)
+    for n ∈ eachindex(tmp)
+        tmp[n] = 0
+        for k ∈ 1 : nbUsers
+            tmp[n] += sig[(n-1)*nbUsers + k] * code[k]
+        end
+        tmp[n] = tmp[n] / nbUsers
+    end
+end
+
+
+# Dispatch with StrucCDMA 
+cdmaSigDecode(sig::AbstractVector,cdma::StrucCDMA) = cdmaSigDecode(sig,cdma.nbUsers,cdma.code,cdma.userMask)