diff --git a/Classes/AWG.m b/Classes/AWG.m index 4c453e3..0ed6ecb 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 = 92e9; + obj.f_cutoff = 32e9; 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 @@ -84,11 +84,11 @@ classdef AWG % Detailed explanation goes here arguments(Input) obj - data_in Informationsignal + data_in end arguments(Output) - elec_out Electricalsignal + elec_out end obj.signal_length = length(data_in); diff --git a/Classes/EQ.m b/Classes/EQ.m new file mode 100644 index 0000000..b52f3f1 --- /dev/null +++ b/Classes/EQ.m @@ -0,0 +1,535 @@ +classdef EQ + %EQ Summary of this class goes here + % Detailed explanation goes here + + properties + Ne %Number of feed forward coefficients (1st, 2nd and 3rd order) + Nb %Number of decision feedback coefficients (1st, 2nd and 3rd order) + K %Number of samples per symbol + delay %Delay of incoming signal + training_length %Number of training symbols + training_loops %Number of loops through sequence for training mode + ideal_dfe %Error free DFE decisions + + DB_aim %Aim at duobinary output sequence + + M %Order of the PAM constellation (only relevant in case of DB aim) + + FFEmu %mu parameter for FFE part in training mode (0 means normalized LMS) + DFEmu % mu parameter for DFE part in training mode + dd_loops % Number of loops through sequence for DD mode + DDmu % mu parameters for DD mode (individual value for each order) + DCmu % mu parameter for the dc tap + + l1act %Activate/deactive l1 regularization + rho %Parameter for speed of coeff shrinking + epsilon %Reciprocal value of the magnitude of the coeff to converge to zero (1st,2nd,3rd order) + thres %Theshold for neglecting coefficienties (1st,2nd,3rd order) + static_act %Activate/deactive static coefficient reduction + + mode2nd %0: no reduction | 1: polynomial | 2: restricted to interval + len_2nd %length of the interval (only for 2nd order mode = 2) + + mode3rd %0: no reduction | 1: polynomial | 2: restricted to interval + len_3rd %length of the interval (only for 3rd order mode = 3/4) + + + %during simulation + k0 + + end + + methods + function obj = EQ(options) + %EQ Construct an instance of this class + % Detailed explanation goes here + arguments(Input) + options + end + + % alles nochmal mappen + end + + function yout = process(obj,data_in,ref_in) + %METHOD1 Summary of this method goes here + % Detailed explanation goes here + + if obj.DB_aim + ref_DB = zeros(size(ref_in)); + for k = 1:length(ref_in) + if k == 1 + ref_DB(k) = ref_in(k); + else + ref_DB(k) = ref_in(k) + ref_in(k-1); + end + end + ref_in = ref_DB; + end + + ref = [zeros(1,obj.Nb(1)-1) ref_in zeros(1,obj.Nb(1))]; + + if isreal(ref) + cplx = 0; + else + cplx = 1; + end + + if obj.static_act + obj.mode2nd = obj.mode2nd + 1; + obj.mode3rd = obj.mode3rd + 1; + else + obj.mode2nd = 1; + obj.mode3rd = 1; + end + + if obj.mode2nd == 1 + N2 = (obj.Ne(2)*(obj.Ne(2)+1))/2; % Number of coefficients for second order + elseif obj.mode2nd == 2 + N2 = obj.Ne(2); + elseif obj.mode2nd == 3 + N2 = (obj.len_2nd+1)*(2*obj.Ne(2)-obj.len_2nd)/2; + elseif obj.mode2nd == 4 + N2 = (ceil(obj.Ne(2)/2)+1)*(2*obj.Ne(2)-ceil(obj.Ne(2)/2))/2; + end + + if obj.mode3rd == 1 + if cplx + N3 = obj.Ne(3)^2*(obj.Ne(3)+1)/2; + else + N3 = obj.Ne(3)*(obj.Ne(3)+1)*(obj.Ne(3)+2)/6; % Number of coefficients for third order + end + elseif obj.mode3rd == 2 + N3 = obj.Ne(3); + elseif obj.mode3rd == 3 + N3 = obj.Ne(3)^2; + elseif obj.mode3rd == 4 + N3 = round(1/6*(obj.len_3rd+1)*(obj.len_3rd+2)*(3*obj.Ne(3)-2*obj.len_3rd)); + elseif obj.mode3rd == 5 + N3 = 2*obj.Ne(3)*obj.len_3rd-obj.len_3rd*(obj.len_3rd+1)+obj.Ne(3); + end + + Nb2 = (obj.Nb(2)*(obj.Nb(2)+1))/2; + Nb3 = obj.Nb(3)*(obj.Nb(3)+1)*(obj.Nb(3)+2)/6; + + data_in = data_in/sqrt(mean(abs(data_in).^2)); % power normalization of input sequence + + if obj.FFEmu == 0 + norm_fac2 = sqrt(mean(abs(data_in.^2).^2)); % power normalization for second and third order terms + norm_fac3 = sqrt(mean(abs(data_in.^3).^2)); % (not necessary, but seems to be more stable if applied --> same as different mu values for linear and nl terms) + else + norm_fac2 = 1; + norm_fac3 = 1; + end + + norm_fac_DFE2 = sqrt(mean(abs(ref_in.^2).^2)); % same for DFE input (reference) + norm_fac_DFE3 = sqrt(mean(abs(ref_in.^3).^2)); + + data = [zeros(1,floor(obj.Ne(1)/2)) data_in zeros(1,obj.Ne(1))]; + + delta_2 = round((obj.Ne(1)-obj.Ne(2))/2); + delta_3 = round((obj.Ne(1)-obj.Ne(3))/2); + + delta_DFE2 = 1;%round((obj.Nb-obj.Nb(2))/2); + delta_DFE3 = 1;%round((obj.Nb-obj.Nb(3))/2); + + % calculate the indices for the combination of second and third order symbols + % - done in advance because it's the same for each iteration, so time + % can be saved + + [ind_mat_2nd,ind_mat_3rd] = obj.calc_ind(obj.Ne(2),N2,obj.Ne(3),N3,obj.mode2nd,obj.mode3rd,obj.len_2nd,obj.len_3rd,cplx); + + [ind_mat_DFE_2nd,ind_mat_DFE_3rd] = obj.calc_DFE_ind(obj.Nb(2),Nb2,obj.Nb(3),Nb3); + + if obj.l1act + epsilon = diag([ones(1,obj.Ne(1))*obj.epsilon(1) ones(1,N2)*obj.epsilon(2) ones(1,N3)*obj.epsilon(3)]); + end + + obj.k0 = obj.delay; % input delay compared to training sequence + + if 1 % obj.active + %% Calculation of the filter coefficients in training based LMS mode + e = zeros(obj.Ne(1)+N2+N3,1); % initialization of filter coefficients + % e(ceil(obj.Ne(1)/2)) = 1; % set central tap to 1 (better starting point since it's closer to the expected solution) + b = zeros(obj.Nb(1)+Nb2+Nb3,1); + e_dc = mean(data_in); % initilaization of the dc tap with the mean value of the data + % e_save = NaN(361,8.6e5); + % save_ind = 1; + for trainloops = 1:obj.training_loops + m = state.k0+1; % starting symbol index at the delay compared to the training sequence + for n = obj.K*state.k0+1:obj.K:obj.K*obj.training_length + m = m+1; + X_1 = data(obj.Ne(1)+n+(obj.K-1):-1:n+obj.K).'; + [X_2,X_3] = calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx); + + D_1 = ref(obj.Nb(1)-state.k0+m-2:-1:m-state.k0-1).'; + [D_2,D_3] = calc_nl_vecs(D_1,ind_mat_DFE_2nd,ind_mat_DFE_3rd,norm_fac_DFE2,norm_fac_DFE3,delta_DFE2,delta_DFE3,cplx); + + input_vec = [X_1;X_2;X_3]; + reference_vec = [D_1;D_2;D_3]; + error = e_dc + e.'*input_vec - b.'*reference_vec - ref_in(m-state.k0); % error = e_dc + e.'*input_vec - b.'*reference_vec - ref_in(m-state.k0); + + if real(obj.FFEmu) + if obj.l1act + sgn_e = e; + sgn_e(e~=0) = e(e~=0)./abs(e(e~=0)); + e = e - obj.rho*sgn_e./(1+epsilon*abs(e)) - error*input_vec*obj.FFEmu; + else + e = e - error*conj(input_vec)*obj.FFEmu; %e = e - error*conj(input_vec)*obj.FFEmu; %e = e - error*input_vec*obj.FFEmu; + end + else + if obj.l1act + sgn_e = e; + sgn_e(e~=0) = e(e~=0)./abs(e(e~=0)); + e = e - obj.rho*sgn_e./(1+epsilon*abs(e)) - error*input_vec/(input_vec.'*input_vec); + else + e = e - error*input_vec/(input_vec.'*input_vec); + end + end + % e_save(:,save_ind) = e; + % save_ind = save_ind+1; + e_dc = e_dc - obj.dcmu*error; + + if obj.Nb(1) > 0 + b = b + obj.DFEmu*error*reference_vec; % Seems like normalized DFE has worse performance + end + end + end + %% + + % Plot the intermediate coefficients after training mode + state.b = b(1:obj.Nb(1)); + state.b2 = b(obj.Nb(1)+1:obj.Nb(1)+Nb2); + state.b3 = b(obj.Nb(1)+Nb2+1:end); + state.e = e(1:obj.Ne(1)); + state.e2 = e(obj.Ne(1)+1:obj.Ne(1)+N2); + state.e3 = e(obj.Ne(1)+N2+1:end); + + if obj.plottrain + figure(8052) + subplot(2,3,1); stem(abs(state.e),'Markersize',2); + title('FFE coeff linear') + xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) + subplot(2,3,2); stem(state.e2,'Markersize',2); + title('FFE coeff nl 2nd') + xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) + subplot(2,3,3); stem(state.e3,'Markersize',2); + title('FFE coeff nl 3rd') + xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) + subplot(2,3,4);stem(state.b,'Markersize',2); + title('DFE coeff linear') + xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) + subplot(2,3,5);stem(state.b2,'Markersize',2); + title('DFE coeff nl 2nd') + xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) + subplot(2,3,6);stem(state.b3,'Markersize',2); + title('DFE coeff nl 3rd') + xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) + set(gcf,'Position',[200 500 700 400]) + end + + if obj.l1act + neg_lin = find(abs(state.e) < obj.thres(1)); + neg_2nd = find(abs(state.e2) < obj.thres(2)); + neg_3rd = find(abs(state.e3) < obj.thres(3)); + + neg = [neg_lin;neg_2nd+obj.Ne(1);neg_3rd+obj.Ne(1)+N2]; % indices of the neglected coefficients + + rel_lin = find(abs(state.e) >= obj.thres(1)); + rel_2nd = find(abs(state.e2) >= obj.thres(2)); + rel_3rd = find(abs(state.e3) >= obj.thres(3)); + + state.coeff_number = length(rel_lin)+2*length(rel_2nd)+3*length(rel_3rd); + + rel = [rel_lin;rel_2nd+obj.Ne(1);rel_3rd+obj.Ne(1)+N2]; % indices of the relevant coefficients + e(neg) = 0; + + ind_mat_2nd(neg_2nd,:) = []; + + ind_mat_3rd(neg_3rd,:) = []; + end + + %% decision directed mode + if ~obj.DB_aim + constellation_in = unique(ref_in); % getting the symbol constellation from reference data + else + if obj.M == 2 + constellation_in = [-3 -2 -1 0 1 2 3]/sqrt(5)*2; + elseif obj.M == 2.5 + constellation_in = [-5 -4 -3 -2 -1 0 1 2 3 4 5]/sqrt(10)*2; + elseif obj.M == 3 + constellation_in = [-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7]/sqrt(21)*2; + else + constellation_in = unique(ref_in); + end + end + + state.constellation_in = constellation_in; + if obj.l1act + coeff = [e(rel);b]; % combine FFE and DFE coefficient vectors for DD mode + else + coeff = [e;b]; + end + + for dd_loop = 1:obj.dd_loops + m = 0; + output_vec = zeros(1,floor(length(data_in)/obj.K)); % initilaization of the output vector + dd_DFE = zeros(obj.Nb(1),1); + D_2 = zeros(Nb2,1); + D_3 = zeros(Nb3,1); + + if all(obj.mu == obj.mu(1)) + mu_mat = obj.mu(1); + else + if obj.l1act + mu_mat = diag([ones(1,length(rel_lin))*obj.mu(1) ones(1,length(rel_2nd))*obj.mu(2) ones(1,length(rel_3rd))*obj.mu(3) ones(1,obj.Nb)*obj.mu(4)]); + else + mu_mat = diag([ones(1,obj.Ne(1))*obj.mu(1) ones(1,N2)*obj.mu(2) ones(1,N3)*obj.mu(3) ones(1,obj.Nb(1)+Nb2+Nb3)*obj.mu(4)]); + end + end + + if obj.load_decisions + pathn = evalin('base','modeldir'); + temp = load([pathn, 'MLSE_out', '.mat']) ; + eval(['dd_out_vals = temp.', 'a', ';']) ; + dd_out = zeros(size(data_in)); + dd_out(1:2:length(data_in)) = dd_out_vals; + else + dd_out = zeros(size(data_in)); + end + + for k = 1:obj.K:length(data_in) + m=m+1; % Symbol index + X_1 = data(obj.Ne(1)+k-1:-1:k).'; + [X_2,X_3] = calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx); + + if obj.l1act + input_vec = [X_1(rel_lin);X_2;X_3;-dd_DFE;-D_2;-D_3]; + else + input_vec = [X_1;X_2;X_3;-dd_DFE;-D_2;-D_3]; + end + + output_vec(m) = e_dc + input_vec.'*coeff; + + if ~obj.load_decisions + [~,dd_idx] = min(abs(output_vec(m) - constellation_in)); % decision for closest constellation point + dd_out(k) = constellation_in(dd_idx); + end + + if obj.Nb(1) > 0 + dd_DFE(2:end) = dd_DFE(1:end-1); + dd_DFE(1) = dd_out(k); + + if obj.error_free && m > state.k0 + dd_DFE(1) = ref_in(m-state.k0); + end + [D_2,D_3] = calc_nl_vecs(dd_DFE,ind_mat_DFE_2nd,ind_mat_DFE_3rd,norm_fac_DFE2,norm_fac_DFE3,delta_DFE2,delta_DFE3,cplx); + end + % if dd_loop ~= 21 + error = output_vec(m) - dd_out(k); + % else + % error = 0; + % end + + coeff = coeff - mu_mat*error*conj(input_vec); + % e_save(:,save_ind) = coeff; + % save_ind = save_ind+1; + + if mu_mat ~= 0 + e_dc = e_dc - obj.dcmu*error; + end + + end + end + + end + + end + + + function [X_2,X_3] = calc_nl_vecs(X_1,ind_mat_2,ind_mat_3,norm_fac2,norm_fac3,delta_2,delta_3,cplx) + % calculation of the vectors containing all combinations of input symbols + % of second and third order based on the linear symbols + + if ind_mat_2(1) > 0 + input_vec_se = X_1(delta_2:end)/norm_fac2;%(K*(k0-1):end) + X_2 = input_vec_se(ind_mat_2(:,1)).*input_vec_se(ind_mat_2(:,2)); + else + X_2 = []; + end + + if ind_mat_3(1) > 0 + if cplx + input_vec_th = X_1(delta_3:end)/norm_fac3; + X_3 = input_vec_th(ind_mat_3(:,1)).*input_vec_th(ind_mat_3(:,2)).*conj(input_vec_th(ind_mat_3(:,3))); + else + input_vec_th = X_1(delta_3:end)/norm_fac3; + X_3 = input_vec_th(ind_mat_3(:,1)).*input_vec_th(ind_mat_3(:,2)).*input_vec_th(ind_mat_3(:,3)); + end + else + X_3 = []; + end + end + + function [ind_mat_2nd,ind_mat_3rd] = calc_ind(Ne2,N2,Ne3,N3,mode2nd,mode3rd,len_2nd,len_3rd,cplx) + + if Ne2 > 0 + ind_mat_2nd = NaN(N2,2); + count=1; + if mode2nd == 1 + for t = 1:Ne2 + for u = t:Ne2 + ind_mat_2nd(count,:) = [t u]; + count = count + 1 ; + end + end + elseif mode2nd == 2 + for t = 1:Ne2 + ind_mat_2nd(t,:) = [t t]; + end + elseif mode2nd == 3 + for t = 1:Ne2 + for u = t:Ne2 + if u-t<=len_2nd + ind_mat_2nd(count,:) = [t u]; + count = count + 1 ; + end + end + end + elseif mode2nd == 4 + for t = 1:Ne2 + for u = t:Ne2 + if u-t<=ceil(Ne2/2) + ind_mat_2nd(count,:) = [t u]; + count = count + 1 ; + end + end + end + end + + else + ind_mat_2nd = 0; + end + + if Ne3 > 0 + ind_mat_3rd = NaN(N3,3); + count=1; + if mode3rd == 1 + if cplx + for t = 1:Ne3 + for u = t:Ne3 + for v = 1:Ne3 + ind_mat_3rd(count,:) = [t u v]; + count = count + 1 ; + end + end + end + else + for t = 1:Ne3 + for u = t:Ne3 + for v = u:Ne3 + ind_mat_3rd(count,:) = [t u v]; + count = count + 1 ; + end + end + end + end + elseif mode3rd == 2 + for t = 1:Ne3 + ind_mat_3rd(t,:) = [t t t]; + end + elseif mode3rd == 3 + for t = 1:Ne3 + for u = t:Ne3 + ind_mat_3rd(count,:) = [t t u]; + if t ~= u + count = count + 1; + ind_mat_3rd(count,:) = [t u u]; + end + count = count + 1; + end + end + elseif mode3rd == 4 + for t = 1:Ne3 + for u = t:Ne3 + for v = u:Ne3 + if u-t<=len_3rd && v-t<=len_3rd + ind_mat_3rd(count,:) = [t u v]; + count = count + 1 ; + end + end + end + end + elseif mode3rd == 5 + for t = 1:Ne3 + for u = t:Ne3 + if u-t<=len_3rd + ind_mat_3rd(count,:) = [t t u]; + if t ~= u + count = count + 1; + ind_mat_3rd(count,:) = [t u u]; + end + count = count + 1; + end + end + end + + ind_mat_3rd2 = NaN(N3,3); + count = 1; + % for t = 1:Ne3 + % ind_mat_3rd2(count,:) = [t t t]; + % count = count + 1; + % end + for t = 1:Ne3 + % ind_mat_3rd2(count,:) = [t t t]; + % count = count + 1; + for u = t:min(Ne3,t+len_3rd) + for v = unique([t u]) + ind_mat_3rd2(count,:) = [t v u]; + count = count + 1; + % ind_mat_3rd2(count,:) = [t u u]; + % count = count + 1; + end + end + end + end + else + ind_mat_3rd = 0; + end + + end + + function [ind_mat_2nd,ind_mat_3rd] = calc_DFE_ind(Ne2,N2,Ne3,N3) + + if Ne2 > 0 + ind_mat_2nd = NaN(N2,2); + count=1; + for t = 1:Ne2 + for u = t:Ne2 + ind_mat_2nd(count,:) = [t u]; + count = count + 1 ; + end + end + else + ind_mat_2nd = 0; + end + + if Ne3 > 0 + ind_mat_3rd = NaN(N3,3); + count=1; + for t = 1:Ne3 + for u = t:Ne3 + for v = u:Ne3 + ind_mat_3rd(count,:) = [t u v]; + count = count + 1 ; + end + end + end + else + ind_mat_3rd = 0; + end + + end + + end +end + diff --git a/comm_tb.m b/comm_tb.m index 47ae742..ce4bd13 100644 --- a/comm_tb.m +++ b/comm_tb.m @@ -17,7 +17,6 @@ fadc = 256e9; % Simulation frequency in "analog domain" fsimu = kover * fdac ; - %SIMULATE for i = 1:log2(M) @@ -32,7 +31,7 @@ Bits.signal = applyPulseShaping(Bits.signal,fsym,fdac); awg = AWG('preset','M8196A','fdac',fdac,'kover',kover,'lpf_active',1); -awgSignal = awg.process_channel(Bits); +awgSignal = awg.process_channel(Bits.signal); fil = Filter('filtdegree',4,"f_cutoff",50e9,"fsamp",fdac,"filterType",filtertypes.bessel_bilin); filtered = fil.process(awgSignal); @@ -62,6 +61,13 @@ scp = Scope("fsimu",fdac*kover,"fadc",fadc,... scope_out = scp.process(phdiod_out2); +resample_out = resample(scope_out,fsym,fadc); + + + + + + figure; hold on plot( nmlze(resample(shapedData,fadc,fdac)),'DisplayName','shapedData Data'); @@ -69,7 +75,7 @@ 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;