diff --git a/Classes/AWG.m b/Classes/AWG.m index 94c7af9..32e5374 100644 --- a/Classes/AWG.m +++ b/Classes/AWG.m @@ -57,7 +57,7 @@ classdef AWG obj.dac_max = 0.5; obj.dac_min = -.5; obj.lowpass = 1; %LP - obj.f_cutoff = 32e9; + obj.f_cutoff = 92e9; elseif options.preset == "M8199B" %https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf end @@ -112,8 +112,10 @@ classdef AWG % 4. Apply LPF on the signal if obj.lpf_active - obj.H_lpf = obj.buildFilter(1); - data_out = obj.lpf(data_out) ; + lpf = Filter('filtdegree',4,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",filtertypes.bessel_inp); + data_out = lpf.process(data_out); +% obj.H_lpf = obj.buildFilter(1); +% data_out = obj.lpf(data_out) ; end end diff --git a/Classes/Filter.m b/Classes/Filter.m index f5e82af..1ea009f 100644 --- a/Classes/Filter.m +++ b/Classes/Filter.m @@ -12,7 +12,7 @@ classdef Filter filtdegree passband_ripple stopband_ripple - fdac + fsamp end methods @@ -22,7 +22,7 @@ classdef Filter arguments options.filterType = 1; options.f_cutoff = 0; - options.fdac = 0; + options.fsamp = 0; options.filtdegree = 3; options.passband_ripple = 0.5; options.stopband_ripple = 0.5; @@ -31,12 +31,12 @@ classdef Filter obj.filterType = options.filterType; - obj.f_cutoff = options.f_cutoff; obj.filtdegree = options.filtdegree; obj.passband_ripple = options.passband_ripple; obj.stopband_ripple = options.stopband_ripple; - obj.fdac = options.fdac; + obj.fsamp = options.fsamp; + end @@ -67,49 +67,49 @@ classdef Filter % Bessel filter, impulse invariant transformed [B, A] = besself(obj.filtdegree, 2*pi*obj.f_cutoff); - [B ,A] = impinvar(B,A,obj.fdac); + [B ,A] = impinvar(B,A,obj.fsamp); case 2 % Bessel filter, impulse bilinear transformed [Z, P, K] = besself(obj.filtdegree, 2*pi*obj.f_cutoff); - [Z ,P, K] = bilinear(Z,P,K,obj.fdac); + [Z ,P, K] = bilinear(Z,P,K,obj.fsamp); [B ,A] = zp2tf(Z ,P ,K); case 3 % Butterworth filter if obj.lowpass == 1 %lowpass - [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fdac/2),'low'); + [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fsamp/2),'low'); else % highpass - [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fdac/2),'high'); + [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fsamp/2),'high'); end case 4 % Chebyshev 1 filter - [B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fdac/2)); + [B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fsamp/2)); case 5 % Chebyshev 2 filter - [B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fdac/2)); + [B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fsamp/2)); case 6 % Elliptic filter - [B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fdac/2)); + [B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fsamp/2)); case 7 % Hamming filter g=(obj.filtdegree-1)/2; - wc=obj.f_cutoff/(obj.fdac/2); + wc=obj.f_cutoff/(obj.fsamp/2); B = wc*sinc(wc*(-g:g)).*hamming(obj.filtdegree)'; A=1; @@ -117,7 +117,7 @@ classdef Filter % Raised Cosine filter - B = firrcos(obj.filtdegree,obj.f_cutoff,para.df,obj.fdac); + B = firrcos(obj.filtdegree,obj.f_cutoff,para.df,obj.fsamp); A=1; case 9 @@ -125,7 +125,7 @@ classdef Filter % Sinc filter g=(obj.filtdegree-1)/2; - wc=obj.f_cutoff/(obj.fdac/2); + wc=obj.f_cutoff/(obj.fsamp/2); B = wc*sinc(wc*(-g:g)); A=1; diff --git a/Classes/Photodiode.m b/Classes/Photodiode.m index 1378061..9d69308 100644 --- a/Classes/Photodiode.m +++ b/Classes/Photodiode.m @@ -37,10 +37,10 @@ classdef Photodiode R = 50; %resistance of phdiode (50ohm is typical value) % Magnitude squared detection - yout = sum( abs(xin) .^2*obj.responsivity,1) ; + yout = sum( abs(xin) .^2*obj.responsivity, 2 ) ; % Shot Noise - shot_noise = sqrt(k * obj.fsimu .* yout) .* randn(1,size(yout,1)); + shot_noise = sqrt(k * obj.fsimu .* yout) .* randn(size(yout,1),1); yout = yout + shot_noise; @@ -50,7 +50,7 @@ classdef Photodiode Bw = obj.fsimu; %is this correct? shouldnt it be the bandwidth of the actual component? e.g. 70GHz? therm_noise_pow = therm_current_psd * 2 * Bw; %squared - therm_noise = sqrt(therm_noise_pow) .* randn(1,size(yout,1)); + therm_noise = sqrt(therm_noise_pow) .* randn(size(yout,1),1); yout = yout + therm_noise; diff --git a/Classes/Scope.m b/Classes/Scope.m new file mode 100644 index 0000000..c3506c3 --- /dev/null +++ b/Classes/Scope.m @@ -0,0 +1,193 @@ +classdef Scope + %SCOPE Summary of this class goes here + % Detailed explanation goes here + + properties + fsimu + fadc + adcresolution + quantbuffer + + rand_samplingdelay %on or off + samplingdelay %specifiy a sampling delay + + freq_offset %offset of the sampler + samp_jitter %include jitter + + fixed_delay %fix the delay of the filter or use minimal delay for kausal system + delay %specify a fixed delay of the filter + + filtertype + lpf_bw + + %during construction + + %during process + Nout + end + + methods + function obj = Scope(options) + %SCOPE Construct an instance of this class + % Detailed explanation goes here + arguments + options.fsimu + options.fadc + options.adcresolution + options.quantbuffer + + options.rand_samplingdelay = 0; + options.samplingdelay = 0; + + options.freq_offset = 0; + options.samp_jitter = 0; + + options.fixed_delay = 0; + options.delay = 0; + + options.filtertype = filtertypes.bessel_bilin; + options.lpf_bw = 120e9; + end + + obj.fsimu = options.fsimu; + obj.fadc = options.fadc; + obj.adcresolution = options.adcresolution; + obj.quantbuffer = options.quantbuffer; + + obj.rand_samplingdelay = options.rand_samplingdelay; % use a randomized sample delay INSTEAD of samplingdelay + obj.samplingdelay = options.samplingdelay; %specifiy a sampling delay + + obj.freq_offset = options.freq_offset; %offset of the sampler + obj.samp_jitter = options.samp_jitter; %include jitter in [s] + + obj.fixed_delay = options.fixed_delay; %fix the delay of the filter or use minimal delay for kausal system + obj.delay = options.delay; %specify a fixed delay of the filter + + obj.filtertype = options.filtertype; + obj.lpf_bw = options.lpf_bw; + + end + + function yout = process(obj,xin) + %METHOD1 Summary of this method goes here + % Detailed explanation goes here + + % apply LPF + lpf = Filter('filtdegree',4,"f_cutoff",obj.lpf_bw,"fsamp",obj.fsimu,"filterType",obj.filtertype); + yout = lpf.process(xin); + + % sample signal + % TODO: implement and test the delays. Also look for delay of + % lpf filter + yout = obj.sampleSignal(yout); + + % quantize signal + yout = obj.quantize(yout); + + + end + + function jitvec = buildSamplingJitterVector(obj) + % Generation of sampling jitter vector for real part of signal + jitvec = obj.samp_jitter*randn(obj.Nout,1) ; + + % Avoid exceding the simulation block edges relentlessly! + jitvec(1) = abs(jitvec(1)) ; + jitvec(end) = -abs(jitvec(end)) ; + + end + + function sampout = sampleSignal(obj,xin) + + % Initialization sampler and output signal in case of 4real + + Tout = 1 / obj.fadc; % output signal sample period [sec] + + Tsig = length(xin) / obj.fsimu ; % simulation block period [sec] + + if obj.rand_samplingdelay == 1 + obj.samplingdelay = Tout*rand ; + end + + %get number of delayed samples + Ndelay = floor(obj.samplingdelay*obj.fsimu) ; + + %get fractional sample delay + fracdelay = obj.samplingdelay - Ndelay/obj.fsimu ; + + %calc number of samples after down-sampling no offset + obj.Nout = ceil((Tsig - fracdelay)*obj.fadc) ; + + %build a jitter vector + jitvec = obj.buildSamplingJitterVector(); + + % Create the index vector of the sampling instances for the real and imaginary part + inx_re = ( (0: Tout : (obj.Nout-1) / obj.fadc)' + jitvec ) * obj.fsimu + 1 ; + + % Create the vector of the relative location of each sample between + % to neighboring "Analog" samples for the real and imaginary part + val_re = mod(inx_re,1) ; + + % Shift the signal + xin = circshift(xin, [0 - Ndelay]) ; + + % Sample the signal using linear interpolation + sampout = xin(floor(inx_re)) .* (1 - val_re) + xin(ceil(inx_re)) .* val_re ; + + end + + function quantout = quantize(obj,xin) + + quant_steps = round(2^obj.adcresolution) + rem(round(2^obj.adcresolution),2); + + if quant_steps > 0 + + index = round(0.1*length(xin)); %sioe 2022: short blocksize resulted in error -> now 10%-90% of length are evaluated + + max_re = max(real(xin(index:end-index)))*(1+obj.quantbuffer/2); % adaptations in the indices + min_re = min(real(xin(index:end-index)))*(1-obj.quantbuffer/2); %(1+para.max/2); + + range_re = max_re-min_re; + + max_im = max(imag(xin))*(1+obj.quantbuffer/2); + min_im = min(imag(xin))*(1+obj.quantbuffer/2); + + range_im = max_im-min_im; + + if isreal(xin) + % shift signal and clip to quantizer intervall + xin = min(max(xin-min_re,0),range_re); + + % quantize signal + xin = round((quant_steps-1)/range_re*xin); + + % scale and shift back + quantout = xin*range_re/(quant_steps-1) +min_re; + + else + + % shift signal and clip to quantizer intervall + xin_re = min(max(real(xin)-min_re,0),range_re); + xin_im = min(max(imag(xin)-min_im,0),range_im); + + % quantize signal + + % quantize signal and shift back signal + xin_re = round((quant_steps-1)/range_re*xin_re); + xin_im = round((quant_steps-1)/range_im*xin_im); + + % scale and shift back + quant_out_re = xin_re*range_re/(quant_steps-1) + min_re; + quant_out_im = xin_im*range_im/(quant_steps-1) + min_im; + quantout = quant_out_re + 1i * quant_out_im; + + end + else + quantout = xin; + end + + end + + end +end + diff --git a/Classes/Signal.m b/Classes/Signal.m index 3c6c202..345897f 100644 --- a/Classes/Signal.m +++ b/Classes/Signal.m @@ -9,13 +9,14 @@ classdef Signal end methods - function obj = Signal(signal, fsym, fsimu) + function obj = Signal(signal, fsym) %SIGNAL Construct an instance of this class % Detailed explanation goes here obj.signal = signal; obj.fsym = fsym; - obj.fsimu = fsimu; + + obj.fsimu = fsym; end @@ -39,8 +40,8 @@ classdef Signal end end - function outputArg = method1(obj,inputArg) - + function normalizedsignal = signal_(obj) + normalizedsignal = obj.signal / norm(obj.signal); end end diff --git a/comm_tb.m b/comm_tb.m index dd876b6..0596170 100644 --- a/comm_tb.m +++ b/comm_tb.m @@ -1,5 +1,5 @@ -O = 10; %order of prbs -N = 64; %length of prbs +O = 14; %order of prbs +N = 2^(O-1); %length of prbs [~,seed] = prbs(O,1); %initialize first seed of prbs % Modulation @@ -12,6 +12,8 @@ fsym = 56e9; fdac = 120e9; % Simulation oversampling rate "k"; kover = 16; +% ADC Rate +fadc = 256e9; % Simulation frequency in "analog domain" fsimu = kover * fdac ; @@ -22,44 +24,73 @@ for i = 1:log2(M) [bitpattern(:,i),seed] = prbs(O,N,seed); end +S = Signal(bitpattern,fsym); + pamData = pam_mapping(bitpattern,M,0); shapedData = applyPulseShaping(pamData,fsym,fdac); - awg = AWG('preset','M8196A','fdac',fdac,'kover',kover,'lpf_active',1); awgSignal = awg.process_channel(shapedData); -fil = Filter('filtdegree',4,"f_cutoff",50e9,"fdac",fdac,"filterType",filtertypes.bessel_bilin); +fil = Filter('filtdegree',4,"f_cutoff",50e9,"fsamp",fdac,"filterType",filtertypes.bessel_bilin); filtered = fil.process(awgSignal); u_pi = 3.5; vbias = (0.5*u_pi)-u_pi; -eml = EML("mode",emlmodes.im_cosinus,"power",10,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1000000); -laserfield = eml.process(filtered); +extmodlaser = EML("mode",emlmodes.im_cosinus,"power",10,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1000000); +laserfield = extmodlaser.process(filtered); att = Amplifier("amplification_db",10,"amp_mode","gain","type","ideal","saturation_mode",0,'saturation_power',10); att_out = att.process(laserfield); -fib = Fiber("fsimu",fdac*kover,"fiber_length",20,"alpha",0.2,"D",17,"lambda0",1550); +fib = Fiber("fsimu",fdac*kover,"fiber_length",0,"alpha",0.2,"D",17,"lambda0",1550); fib_out = fib.process(att_out); phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); phdiod_out = phdiode.process(fib_out); +fil = Filter('filtdegree',4,"f_cutoff",70e9,"fsamp",fdac,"filterType",filtertypes.bessel_bilin); +phdiod_out2 = fil.process(phdiod_out); +scp = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"lpf_bw",120e9,"filtertype",filtertypes.bessel_inp,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",6,"quantbuffer",0.1); + +scope_out = scp.process(phdiod_out2); + figure; -plot(shapedData); hold on -plot(awgSignal,'DisplayName','skew 0'); -plot(filtered,'DisplayName','filtered'); -plot(abs(laserfield),'DisplayName','laser'); -plot(abs(att_out),'DisplayName','att_out'); -plot(abs(fib_out),'DisplayName','att_out'); - +plot( nmlze(resample(shapedData,fadc,fdac)),'DisplayName','shapedData Data'); +plot( nmlze(resample(awgSignal,fadc,fdac*kover)),'DisplayName','AWG Data'); +plot( nmlze(scope_out),'DisplayName','scope out out'); hold off +resample_out = resample(scope_out,fsym,fadc); + + +% figure; +% hold on +% plot(awgSignal,'DisplayName','skew 0'); +% plot(filtered,'DisplayName','filtered'); +% plot(abs(laserfield),'DisplayName','laser'); +% plot(abs(att_out),'DisplayName','att_out'); +% plot(abs(fib_out),'DisplayName','fiber_out'); +% plot(abs(phdiod_out),'DisplayName','photo diode out'); +% hold off + +% figure; +% hold on +% plot(nmlze(pamData),'DisplayName','PAM Data'); +% plot(nmlze(resample_out),'DisplayName','system out'); +% hold off + + +function y = nmlze(x) + y = (x-min(x)) / max((x-min(x))); +end function pam_sig = pam_mapping(bitpattern, M, unipolar)