Mainly AWG updates

Classes/01_transmit/AWG.m

@@ -1,4 +1,4 @@
-classdef AWG
+classdef AWG < handle
     %AWG Summary of this class goes here
     %   Detailed explanation goes here
 
@@ -63,6 +63,7 @@ classdef AWG
             len_in = length(signalclass_in.signal);
 
             if signalclass_in.fs ~= obj.fdac
+                k = obj.kover*obj.fdac / signalclass_in.fs;
                 signalclass_in = signalclass_in.resample("fs_in",signalclass_in.fs,"fs_out",obj.fdac);
             end
                        
@@ -70,11 +71,14 @@ classdef AWG
             signalclass_in.signal = obj.process_(signalclass_in.signal);
             
             % cast the inform. signal to electrical signal
+            if ~isa(signalclass_in,'Electricalsignal')
                 signalclass_in = Electricalsignal(signalclass_in,"fs",obj.fdac*obj.kover,"logbook",signalclass_in.logbook);
+            end
 
             % normalize to 0dBm before applying the lowpass
             %signalclass_in = signalclass_in.normalize("mode","milliwatt");
-            signalclass_in = signalclass_in.setPower(6,"dBm");
+            signalclass_in = signalclass_in.setPower(13,"dBm");
+
 
             % 4. Apply LPF on the signal
             if obj.lpf_active
@@ -89,7 +93,7 @@ classdef AWG
                 end
             end
             
-            signalclass_in.power
+            disp(['AWG output power: ',num2str(signalclass_in.power),' dBm']);
             
             % append to logbook
             current_class = class(obj);
@@ -123,12 +127,15 @@ classdef AWG
 
             if obj.normalize2dac
                 % 0a Normalize the signal to full scale DAC range
-                data_in = data_in - min(data_in);
+               
-                data_in = data_in/(max(data_in)-min(data_in));
+                data_in = data_in - min(data_in); %set "foot" to zero
-                data_in = data_in * (obj.dac_max-obj.dac_min);
+                data_in = data_in / (max(data_in)-min(data_in)); %scale between 0 and 1
-                data_in = data_in + obj.dac_min;
+                data_in = data_in * (obj.dac_max-obj.dac_min); %scale to desired total range of DAC
+                data_in = data_in + obj.dac_min; %set "foot" to desired value
             end
 
+            data_in = data_in - mean(data_in);
+
             % 1. Quantize the signal - Full Scale is between obj.dac_min
             % and dac_max. If signal is smaller in between, you won't use
             % the full bit-resolution.
@@ -138,6 +145,7 @@ classdef AWG
                 elec_out = data_in;
             end
              
+
             % 2. Sample and hold + repeat (data_out: 1xsignal length)
             if obj.upsampling_method == 1
                 % just use matlab function
@@ -150,6 +158,7 @@ classdef AWG
                 error('chosen upsampling method not implemented?');
             end
           
+
             % 3. Add skew (not implemented so far)
             if obj.skew_active
                 elec_out = obj.skew(elec_out);
@@ -182,7 +191,7 @@ classdef AWG
                     1i*min(max(imag(x_in)-obj.dac_min,0),obj.dac_max-obj.dac_min);
 
                 % quantize signal and shift back signal
-                x_in = round((steps-1)/(obj.max-obj.dac_min)*x_in);
+                x_in = round((steps-1)/(obj.dac_max-obj.dac_min)*x_in);
 
                 % scale and shift back
                 quant_out = x_in*(obj.dac_max-obj.dac_min)/(steps-1) +(1+1i)*obj.dac_min;

Classes/02_etc/Filter.m

@@ -249,8 +249,12 @@ classdef Filter < handle
         end
 
         function showHere(obj)
-            obj.signal_length = 512;
+            if isempty(obj.H)
+                obj.signal_length = 4096;
                 [H_,~] = obj.buildFilter();
+            else
+                H_ = obj.H;
+            end
             fs = obj.fs;
             frex = linspace(-fs/2,fs/2,numel(H_)).*1e-9;
 
@@ -265,6 +269,12 @@ classdef Filter < handle
             [~, index] = min(abs(filt - (-3)));
             threedB = frex(index);
 
+            [~, index] = min(abs(filt - (-6)));
+            sixdB = frex(index);
+
+            [~, index] = min(abs(filt - (-10)));
+            ninedB = frex(index);
+
             cur_module_name = '';
 
             filter_desc = char([num2str(obj.filtdegree),'th-order ',char(obj.filterType),'; f cut:',num2str(fc.*1e-9),' GHz']);
@@ -273,10 +283,17 @@ classdef Filter < handle
             legend Interpreter none
             title('Magnitude')
             hold on
-            p = plot(frex,filt,'LineWidth',1,'LineStyle','-','DisplayName',[cur_module_name,': ',filter_desc]);
+            p = plot(frex,filt,'LineWidth',1,'LineStyle','-','DisplayName',[cur_module_name,': ',filter_desc, '; 3/6/10 dB: ', num2str(threedB),'/ ',num2str(sixdB),'/ ',num2str(ninedB)]);
+
             xline([-threedB,threedB],'LineStyle','-','LineWidth',1,'HandleVisibility','off','Color',p.Color);
+            xline([-sixdB,sixdB],'LineStyle','-','LineWidth',1,'HandleVisibility','off','Color',p.Color);
+            xline([-ninedB,ninedB],'LineStyle','-','LineWidth',1,'HandleVisibility','off','Color',p.Color);
+
             xline(fcut_,'LineStyle',':','LineWidth',1,'HandleVisibility','off','Color',p.Color);
-            yline(-3,'LineStyle',':','LineWidth',1,'HandleVisibility','off');
+            yline([-3, -6, -9],'LineStyle',':','LineWidth',1,'HandleVisibility','off');
+            
+            xlim([0 fc.*2].*1e-9);
+            ylim([ninedB-6, 2]);
             legend
 
         end

Functions/plot_eye.m

@@ -19,6 +19,9 @@ function plot_eye(signal, fs, fsym)
     maxA = max(eye_mat(:))*1.1;
     minA = min(eye_mat(:))*1.1;
 
+    maxA = max(x);
+    minA =min(x);
+
     difference= maxA-minA;
 
     data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1;
@@ -28,7 +31,7 @@ function plot_eye(signal, fs, fsym)
         hist_data(:,n)=hist_data(:,n)+nn.';
     end
 
-    figure;
+
     image(hist_data);
     colormap(cbrewer2("Reds",92));
 

projects/MPI_April/mpi_simulation_cspr.m

@@ -1,7 +1,7 @@
 %% Settings
 
 clear
-for M=[4,6,8]
+for M=[4]
 
 
     filename = '112G_2';
@@ -11,7 +11,7 @@ for M=[4,6,8]
 
     kover = 8;
     fsym = round(datarate*1e-9 / log2(M))*1e9;
-    fdac = 256e9;
+    fdac = fsym;%256e9;
     fadc = 256e9;
 
     lowpass_cutoff = fsym/2 * 1.1;
@@ -20,21 +20,22 @@ for M=[4,6,8]
     phd_bw = lowpass_cutoff;
     scp_bw = lowpass_cutoff;
 
-    LP_awg      = Filter('filtdegree',4,"f_cutoff",90e9,"fs",fdac*kover,"filterType",filtertypes.butterworth);
+    LP_awg      = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*kover,"filterType",filtertypes.butterworth);
     LP_modulator= Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.butterworth);
     LP_opt      = Filter('filtdegree',3,"f_cutoff",fsym/log2(M).*1.5,"fs",fdac*kover,"filterType",filtertypes.gaussian);
     LP_phd      = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.butterworth);
     LP_scpe     = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth);
 
-
     % 1) PRBS Generation
     O = 18; %order of prbs
     N = 2^(O-1); %length of prbs
     [~,seed] = prbs(O,1); %initialize first seed of prbs
     bitpattern=[];
+
     for i = 1:log2(M)
         [bitpattern(:,i),seed] = prbs(O,N,seed);
     end
+
     if M == 6
         bitpattern = reshape(bitpattern,[],1);
         bitpattern = bitpattern(1:end-mod(length(bitpattern),5));
@@ -62,7 +63,7 @@ for M=[4,6,8]
     vp = [0.25,0.5,0.75,1];
     vb = [1:0.1:1.8];
 
-    rop = -5:5;
+    rop = 0;
     sir = 28;
     laser_linewidth = 10e6;
     pn_key = 9;
@@ -81,10 +82,59 @@ for M=[4,6,8]
                         %digimod_out = digimod_out.normalize("mode","oneone");
 
                         % cnt = cnt+1;
-                        X = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",16,"rrcalpha",0.5).process(digimod_out);
+                        %X = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",16,"rrcalpha",0.05).process(digimod_out);
+                        X = digimod_out;
+
+%                         % precomp 
+                        [b,a] = butter(2,75e9/(fsym/2),"low");
+                        H=freqz(b,a,length(X),fsym,'whole');
+                        max_amp_db = 3;
+                        p = find(abs(H)<10^(-max_amp_db/20));
+                        H(p) = 10^(-max_amp_db/20);
+                        H_maxatt = H;
+                        H_maxatt = max(H,10^(-20/20));
+                        H_inv = 1./H_maxatt;
+                        
+                        freq_vec = linspace(-X.fs/2,X.fs/2,length(X));
+
+                        h = sin(pi*freq_vec/fdac)./(pi*freq_vec/fdac);
+                        
+                        max_amp_db = 40;
+                        max_amp_lin = 10^(-max_amp_db/20);
+
+                        p = find(h<max_amp_lin);
+
+                        h(p)=10^(-max_amp_db/20);
+
+                        h = 1./h;
+
+                        figure(13);
+                        plot(freq_vec.*1e-9,20*log10(fftshift(H_inv)));
+                        hold on
+                        plot(freq_vec.*1e-9,20*log10((h)));
+                        
+
+                         
+                         X.signal = ifft(fftshift(h).'.*fft(X.signal));
+%                         X.signal = ifft(H_inv.*fft(X.signal));
+
+                        figure(12)
+                        spectrum_plot(X.signal,fsym);                 
   
                         % 5) AWG (lowpass, quantization, sample and hold)
-                        X = AWG("fdac",fdac,"dac_min",-1,"dac_max",1,"H_lpf",LP_awg,"kover",kover,"bit_resolution",5,"lpf_active",1).process(X);
+                        X = AWG("fdac",fdac,"dac_min",-1,"dac_max",1,"lpf_active",0,"H_lpf",LP_awg,"kover",kover,"bit_resolution",16,"normalize2dac",1,"upsampling_method","samplehold").process(X);
+                            
+                        spectrum_plot(X.signal,X.fs);
+                           
+
+                        burg_coeff = arburg(X.signal(10000:20000),100);
+
+                            [h,w] = freqz(1,burg_coeff,length(X),"whole",fsym*kover);
+                            h = h/max(abs(h));
+                            hold on
+                            plot(w,20*log10(h),'DisplayName','Burg Coeff');
+
+                        spectrum_plot(X.signal,X.fs);
 
                         % 6) Lowpass behavior before laser
                         X =  LP_modulator.process(X);
@@ -100,11 +150,11 @@ for M=[4,6,8]
 
                         [Opt,extmodlaser] = extmodlaser.process(E);
                         
-                        %                 figure(m)
+                        figure(m)
-                        %                 hold on
+                        hold on
-                        %                 scatter(E.signal(1:100000),(abs(Opt.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF')
+                        scatter(E.signal(1:100000),(abs(Opt.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF')
-                        %                 xlabel('Input in V')
+                        xlabel('Input in V')
-                        %                 ylabel('abs(Output) in mW')
+                        ylabel('abs(Output) in mW')
 
                         % ER = 10*log10(max(abs(Opt.signal).^2)/min(abs(Opt.signal).^2));
 
@@ -131,7 +181,7 @@ for M=[4,6,8]
                         % Fiber
                         Combined_sig    = Fiber("fsimu",Combined_sig.fs,"fiber_length",2,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.08).process(Combined_sig);
 
-                        parfor i = 1:length(rop)
+                        for i = 1:length(rop)
 
                             % Set ROP
                             Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",rop(i)).process(Combined_sig);
@@ -181,16 +231,8 @@ for M=[4,6,8]
                                     "eq_avg_blocklength",0).process(Scpe_sig,digimod_out);
 
                             % Demap
-%                             Rx_Bits = PAMmapper(M,0).demap(EQ_sig);
+
                             Rx_Bits = PAMmapper(M,0).demap(EQ_sym);
-                            %
-                            %                         Rx_symboldecision = PAMmapper(M,0).decide_pamlevel(EQ_sig,"symbol_levels",unique(digimod_out.signal));
-                            %
-                            %                         levels =    PAMmapper(M,0).separate_pamlevels(EQ_sig);
-                            %                         %levels =    PAMmapper(M,0).separate_pamlevels(Rx_sig.resample("fs_in",Rx_sig.fs,"fs_out",fsym));
-                            %
-                            %                         level_avg(s,l,pnk,n,m,i,:) = mean(levels,'omitnan');
-                            %                         level_std(s,l,pnk,n,m,i,:) = std(levels,'omitnan');
 
                             % BER
                             [~,errors_bm,BER(s,l,pnk,n,m,i),errors] = calc_ber(Rx_Bits.signal,bitpattern,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
@@ -198,13 +240,13 @@ for M=[4,6,8]
                             formatted_ber = sprintf('%.1e', BER(s,l,pnk,n,m,i));
                             disp(['SIR: ',num2str(sir(s)),'; Lw:',num2str(laser_linewidth(l)),'; Key:',num2str(pn_key(pnk)),'; Vpeak: ',num2str(vp(n)),'; Vbias',num2str(vbias),'; BER: ',formatted_ber,'; run: ',num2str(cnt),' / 12961']);
 
-                            %                          plot_analysis_window;
+                            plot_analysis_window;
                             %                          drawnow;
 
                         end
                     end
                 end
-                save('pam4_level_comp');
+                % save('pam4_level_comp');
                 disp('saved_run2');
             end
         end

projects/MPI_April/simulate_tx_sig.m

@@ -14,8 +14,6 @@ function simulate_tx_sig(filename)
     
     lp_laser    = Filter('filtdegree',2,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.butterworth);
     
-
-    
     %% PROCESS TX
     
     % 1) PRBS Generation

test/modulator_test.m

@@ -3,7 +3,7 @@
 M = 4;
 laser_linewidth = 0;
 kover = 32;
-fsym = 92e9;%round(datarate*1e-9 / log2(M))*1e9;
+fsym = 170e9;%round(datarate*1e-9 / log2(M))*1e9;
 fdac = 256e9;
 
 % 1) PRBS Generation
@@ -26,41 +26,44 @@ digimod_out = PAMmapper(M,0).map(bits);
 digimod_out.fs = fsym;
 
 % 3) Pulseform Raised Cosine
-X = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",16,"rrcalpha",0.2).process(digimod_out);
+X = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",16,"rrcalpha",0.01).process(digimod_out);
 
 % Implememt Precompensation
 
 % Implement Precoding
 
 % 4) AWG (lowpass, quantization, sample and hold)
-LP_awg = Filter('filtdegree',5,"f_cutoff",64e9,"fs",fdac*kover,"filterType",filtertypes.butterworth);
+
+LP_awg = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*kover,"filterType",filtertypes.gaussian);
+
+
 AWG_=AWG("fdac",fdac,"dac_min",-1,"dac_max",1,"H_lpf",LP_awg,"kover",kover,"bit_resolution",16,"lpf_active",1,"normalize2dac",1,"upsampling_method","samplehold");
 X = AWG_.process(X);
-plot_eye(X.signal,X.fs,fsym);
+
 disp(['El. power: ',num2str(X.power),' dBm (into 50 Ohm)']);
 disp(['El. RMS voltage: ',num2str(sqrt(mean(X.signal.^2))),' V']);
 disp(['max voltage: ',num2str(max(X.signal)),' V']);
 
+
 % 5) Lowpass behavior before laser
 LP_modulator= Filter('filtdegree',4,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.butterworth);
 X =  LP_modulator.process(X);
 
 % 6) Laser; Modulation ->  OPTICAL DOMAIN
 u_pi = 4;
-vbias = -1.8;
+vbias = 2;
 extmodlaser = EML("mode",eml_mode.im_cosinus,"power",0,"fsimu",X.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",5);
 [Opt,extmodlaser] = extmodlaser.process(X);
 
 if 1
-    f = figure(22);
+    f = figure(120);
-    tiledlayout(2,2);
+    f.Name = 'bla';
-    nexttile
+    tiledlayout(2,4);
-    hold on
-    scatter(X.signal(1:100000),(abs(Opt.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF')
-    xlabel('Input in V')
-    ylabel('Output in mW')
 
     nexttile
+    rms_ = rms(X.signal);
+    max_ = max(X.signal);
+    min_ = min(X.signal);
     hold on
     plot(X.signal,'LineWidth',0.1);
     yline([max_, min_],'LineWidth',2,'LineStyle','--');
@@ -68,16 +71,38 @@ if 1
     ylim([-3 3]);
     title(['AWG output: ',num2str(X.power), 'dBm']);
 
+    % Add text boxes for MIN, MAX, and RMS voltage
+    text(0.5, min_-0.3, ['MIN: ', num2str(min_),' V'],'FontSize', 10, 'HorizontalAlignment', 'left');
+    text(0.5, max_+0.3, ['MAX: ', num2str(max_),' V'],'FontSize', 10, 'HorizontalAlignment', 'left');
+    text(0.5, rms_+0.22, ['RMS: ', num2str(rms_),' V'],'FontSize', 10, 'HorizontalAlignment', 'left');
+    text(0.5, -rms_-0.22, ['RMS: ', num2str(rms_),' V'],'FontSize', 10, 'HorizontalAlignment', 'left');
 
     nexttile
-    spectrum_plot(Opt.signal,Opt.fs);
+    plot_eye(X.signal,X.fs,fsym);
-    
+    ylabel('Signal in V')
-    
+
-    rms_ = rms(X.signal);
+    nexttile
-    max_ = max(X.signal);
+    hold on
-    min_ = min(X.signal);
+    v_in_curve = [-u_pi*1.5/2:0.1:u_pi*1.5/2];
+    field=sqrt(10^(extmodlaser.power/10-3));
+    mzm_curve = ((field.*cos(pi/2*(real(v_in_curve)+vbias)/u_pi)).^2)*1e3;
+    scatter(v_in_curve+vbias,mzm_curve,10,'o','filled','DisplayName','Modulator TF complete');
+    scatter(X.signal(1:100000)+vbias,(abs(Opt.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF')
+    scatter(min_+vbias,((field.*cos(pi/2*(real(min_)+vbias)/u_pi)).^2)*1e3,50,'x','LineWidth',2);
+    scatter(max_+vbias,((field.*cos(pi/2*(real(max_)+vbias)/u_pi)).^2)*1e3,50,'x','LineWidth',2);
+    xlim([-u_pi*1.5/2+vbias, u_pi*1.5/2+vbias]);
+    ylim([min(mzm_curve),max(mzm_curve)]);
+    xlabel('Input in V')
+    ylabel('Output in mW')
+    title("MZM input (v) to output (w)");
+
+    nexttile
+    plot_eye(abs(Opt.signal.^2).*1e3 ,Opt.fs,fsym);
+    ylabel('Opt. Signal in mW')
     
 
+    nexttile([1 2])
+    spectrum_plot( Opt.signal,Opt.fs, 'bla');
 end