diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 8224f4b..cb5c411 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -337,15 +337,24 @@ classdef Signal obj options.fignum options.displayname = ""; + options.color = []; + options.normalizeToNyquist = 0; + options.normalizeTo0dB = 0; end % spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname); N = 2^(nextpow2(length(obj.signal))-8); - [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); - normalize = 0; - if normalize + if options.normalizeToNyquist==0 + [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); + w=w.*1e-9; + else + [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,"centered","power","mean"); + end + + + if options.normalizeTo0dB p_lin = p_lin./ max(p_lin); p_dbm = 10*log10(p_lin); %dB to dBm in case of "power" ylab = "normalized to 0 dB"; @@ -357,14 +366,27 @@ classdef Signal figure(options.fignum); % If figure does not exist, create new figure ax = gca; hold on - plot(w.*1e-9,p_dbm,'DisplayName',options.displayname,'LineWidth',1); - xlabel("Frequency in GHz"); - %ylabel("Power/frequency (dB/Hz)"); + + if isempty(options.color) + plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1); + else + plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color); + end + + if options.normalizeToNyquist==0 + xlabel("Frequency in GHz"); + %xlim([-obj.fs/2 obj.fs/2].*1e-9) + edgetick = 2^(nextpow2(obj.fs*1e-9)); + xticks([-edgetick:16:edgetick]); + xlim([100*round( min(w.*1e-9)/100,1)-10,100*round( max(w.*1e-9)/100,1)+10]) + else + xlabel("Normalized Frequency"); + xlim([-pi, pi]); + end + + ylabel("Power/frequency (dB/Hz)"); ylabel(ylab); - xlim([-obj.fs/2 obj.fs/2].*1e-9) - edgetick = 2^(nextpow2(obj.fs*1e-9)); - % xticks([-edgetick:16:edgetick]); - xlim([100*round( min(w.*1e-9)/100,1)-10,100*round( max(w.*1e-9)/100,1)+10]) + ylim([min(floor( min(p_dbm))-3 , ax.YLim(1)), max(ceil( max(p_dbm) )+(3), ax.YLim(2))]); yticks([-200:10:10]); grid on diff --git a/Classes/04_DSP/Equalizer/EQ.m b/Classes/04_DSP/Equalizer/EQ.m index 888f85a..f5bf2a4 100644 --- a/Classes/04_DSP/Equalizer/EQ.m +++ b/Classes/04_DSP/Equalizer/EQ.m @@ -104,7 +104,7 @@ classdef EQ end - function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in) + function [signalclass_out,noi] = process(obj,signalclass_in, reference_signalclass_in) % actual processing of the signal (steps 1. - 3.) [signalclass_in.signal,error_log] = obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); @@ -120,6 +120,9 @@ classdef EQ % write to output signalclass_out = signalclass_in; + noi = signalclass_out; + noi = signalclass_out - reference_signalclass_in; + end function [yout,error_log] = process_(obj,data_in,ref_in) diff --git a/Classes/05_Lab/AwgKeysight.m b/Classes/05_Lab/AwgKeysight.m index 0270ed4..912eb10 100644 --- a/Classes/05_Lab/AwgKeysight.m +++ b/Classes/05_Lab/AwgKeysight.m @@ -215,7 +215,7 @@ classdef AwgKeysight switch obj.model case awg_model.M8196A v = visadev('TCPIP0::zizou.tf.uni-kiel.de::hislip0::INSTR'); - case awg_model.M8196A + case awg_model.M8199B v = visadev('TCPIP0::localhost::hislip0::INSTR'); end @@ -224,10 +224,10 @@ classdef AwgKeysight disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); end - %check if channel config matches - writeline(v,'*opt?'); - aw=readline(v); - assert(obj.numChannels==sscanf(aw, '%f'),['Wrong Channel config. AWG response: ',num2str(sscanf(aw, '%f')), ' Channels']); + % %check if channel config matches (ONLY FOR M8196...) + % writeline(v,'*opt?'); + % aw=readline(v); + % assert(obj.numChannels==sscanf(aw, '%f'),['Wrong Channel config. AWG response: ',num2str(sscanf(aw, '%f')), ' Channels']); switch obj.model case awg_model.M8196A @@ -276,33 +276,27 @@ classdef AwgKeysight case awg_model.M8199B - % 1) TRAC DEl - delete the previous waveform - xfprintf(f, sprintf(':TRAC:DEL ''M2.DataOut1'', %d', segm_num), 1); - xfprintf(f, sprintf(':TRAC:DEL ''M2.DataOut2'', %d', segm_num), 1); + % 1) TRAC DEl - delete the previous ch waveform + obj.writeNcheck(v,sprintf(':TRAC:DEL ''M2.DataOut%d'', %d', chan, segm_num)); % 2) TRAC DEF - define properties of waveform - xfprintf(f, sprintf(':TRAC:DEF ''M2.DataOut1'', %d,%d', segm_num, len)); - xfprintf(f, sprintf(':TRAC:DEF ''M2.DataOut2'', %d,%d', segm_num, len)); + obj.writeNcheck(v,sprintf(':TRAC:DEF ''M2.DataOut%d'', %d,%d', chan, segm_num, segm_len)); + % xfprintf(f, sprintf(':TRAC:DEF ''M2.DataOut1'', %d,%d', segm_num, len)); % 3) TRAC DATA - load the new waveform % Example command: :TRAC:DATA 'M1.DataOut1', 1, 0 , -127,-126,,-1,0,1,,127,-127,,-1,0,1,,127 % means : :TRAC:DATA 'identifier' , , , | - - cmd = sprintf(':TRAC:DATA ''M2.DataOut1'', %d,%d,', segm_num, offset); - binblockwrite(f, data(1+offset:offset+len), dataFormat, cmd); - - cmd = sprintf(':TRAC:DATA ''M2.DataOut2'', %d,%d,', segm_num, offset); - binblockwrite(f, data(1+offset:offset+len), dataFormat, cmd); + offset = 0; + cmd = sprintf(':TRAC:DATA ''M2.DataOut%d'', %d,%d,',chan, segm_num, offset); + binblockwrite(v, signal{chan}, 'int8', cmd); % 4) SET VOLTAGE (not possible to set 0v) - xfprintf(f, sprintf(':VOLT:AMPL ''M2.DataOut1'', %gv', awg_volt(1))); - xfprintf(f, sprintf(':VOLT:AMPL ''M2.DataOut2'', %gv', awg_volt(2))); + obj.writeNcheck(v, sprintf(':VOLT:AMPL ''M2.DataOut%d'', %gv',chan, obj.voltages(chan))); % xfprintf(f, ':VOLT:AMPL ''M2.DataOut2'', 0.1v'); % 5) SET OUTPUTS ON/ OFF - xfprintf(f, ':OUTP ''M2.DataOut1'', ON'); %sioe 11.23: Hier habe ich absichtlich auf OFF gestellt, wenn Kanal nicht bentigt! - xfprintf(f, ':OUTP ''M2.DataOut2'', ON'); + obj.writeNcheck(v, sprintf(':OUTP ''M2.DataOut%d'', ON',chan)); %sioe 11.23: Hier habe ich absichtlich auf OFF gestellt, wenn Kanal nicht bentigt! @@ -319,6 +313,9 @@ classdef AwgKeysight case awg_model.M8199B + % SET Ref clock out (which is conn. to scope) to 100 MHz + obj.writeNcheck(v, ':OUTPut:FREQuency ''M1.RefClkOut'',100000000'); + % :OUTP:GLOB 'M1.System', OFF obj.writeNcheck(v, ':OUTP:GLOB ''M1.System'', ON'); obj.writeNcheck(v, ':OUTP:GLOB ''M2.System'', ON'); diff --git a/Classes/05_Lab/DC_supply.m b/Classes/05_Lab/DC_supply.m index a04cd3e..5c17d68 100644 --- a/Classes/05_Lab/DC_supply.m +++ b/Classes/05_Lab/DC_supply.m @@ -1,4 +1,4 @@ -classdef DC_supply +classdef DC_supply < handle % Simple control of DC supply via fixed GBIB address. % Agilent E3646A % No OVP or OCP implementation @@ -8,6 +8,10 @@ classdef DC_supply properties(Access=public) active voltage + + + meas_voltage + meas_current end methods (Access=public) @@ -33,6 +37,10 @@ classdef DC_supply end + function v = connectDevice(obj) + v = visadev("GPIB1::19::INSTR"); + end + function success = set(obj,options) @@ -44,7 +52,14 @@ classdef DC_supply success = [0,0]; try - %connect to device + % %connect to device + % if exist('v','var') + % %v = visadev("GPIB1::19::INSTR"); + % v = obj.connectDevice(); + % else + % + % end + v = visadev("GPIB1::19::INSTR"); debug = 0; @@ -121,5 +136,55 @@ classdef DC_supply end end + + function [voltage,current] = readVals(obj) + + try + %connect to device + v = visadev("GPIB1::19::INSTR"); + + debug = 0; + if debug + disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); + end + + cmd = 'INST:SEL?'; + writeline(v, cmd); + prev_selected_channel = readline(v); + + for ch = 1:2 + + % choose channel + cmd = ['INST:SEL OUT',num2str(ch)]; + writeline(v, cmd); + + % get current voltage level + cmd = 'VOLT?'; + writeline(v, cmd); + voltage(ch) = str2num(readline(v)); + + % get current voltage level + cmd = 'MEAS:CURR?'; + writeline(v, cmd); + current(ch) = str2num(readline(v)); + + + end + + % set channel back to prev. selected one + cmd = ['INST:SEL ',char(strtrim(prev_selected_channel))]; + writeline(v, cmd); + + obj.meas_voltage = voltage; + obj.meas_current = current; + + %disconnect + delete(v); + + catch + + end + end + end end diff --git a/projects/Lab_2024/lab_modulator_tf_sweep.m b/projects/Lab_2024/lab_modulator_tf_sweep.m index 6c41e4a..2a04cec 100644 --- a/projects/Lab_2024/lab_modulator_tf_sweep.m +++ b/projects/Lab_2024/lab_modulator_tf_sweep.m @@ -1,68 +1,136 @@ +measurement_name = 'MZM_SD40'; + % Initialize a structure to hold parameters params = struct; % Define the bias voltage range from 1.2V to 2.8V with 0.01V increments -params.v_bias = 1.2:0.01:2.8; +params.v_bias = 0:0.1:6; % Create a DataStorage object with the defined parameters wh = DataStorage(params); % Add a storage field for output power measurements wh.addStorage("p_out"); +wh.addStorage("i_draw"); +wh.addStorage("v_set"); % Display the total number of measurement loops to be executed disp(['Start Measurement of ', num2str(prod(wh.dim)), ' loops...']); + % Initialize the Optical Attenuator (VOA) with specified settings voa = OptAtten(... - "active", [1, 0, 0, 0], ... % Activate only the first channel + "active", [0, 1, 0, 0], ... % Activate only the first channel "value", [0, 0, 0, 0], ... % Set attenuation values to 0 dB "wavelength", [1310, 1310, 1310, 1310]); % Set the wavelength for each channel +% Set the VOA active channels and attenuation values +voa.set('active', [0, 1, 0, 0], 'value', [0, 0, 0, 0]); + % Initialize the DC Power Supply with specified settings -dcs = DC_supply(... - "active", [1, 1], ... % Activate the first two channels - "voltage", [v_bias, 9]); % Set initial voltages for channels +v_bias_init = wh.parameter.v_bias.values(1); + +dcs = DC_supply(... + "active", [1, 0], ... % Activate the first two channels + "voltage", [v_bias_init, 0]); % Set initial voltages for channels -% Set the VOA active channels and attenuation values -voa.set('active', [1, 0, 0, 0], 'value', [0, 0, 0, 0]); % Loop over each bias voltage value to perform measurements for v_bias = wh.parameter.v_bias.values - try - % Try to retrieve existing output power measurement to avoid repetition - p_out = wh2.getStoValue('p_out', v_bias); - wh.addValueToStorage(p_out, 'p_out', v_bias); - catch - % If no existing measurement, proceed with the measurement - - %%%%% SET Voltages %%%%%% - % Update the DC supply voltage for the current bias voltage - dcs.set("voltage", [v_bias, 9]); - - %%%%% Measure VOA %%%%%% - % Read the current values from the VOA - voa.readvals(); - - % Store the measured output power in the DataStorage object - wh.addValueToStorage(voa.power_state(1), 'p_out', v_bias); - end + + %%%%% SET Voltages %%%%%% + % Update the DC supply voltage for the current bias voltage + dcs.set("voltage", [v_bias, 0]); + + %%%%% Measure Voltages %%%%%% + [v,c] = dcs.readVals(); + % assert(v(1)==v_bias,['Desired voltage is ', num2str(v_bias),'V but currently is ', num2str(v(1)), ' V ']); + + %%%%% Measure VOA %%%%%% + % Read the current values from the VOA + voa.readvals(); + + % Store the measured output power in the DataStorage object + wh.addValueToStorage(voa.power_state(2), 'p_out', v_bias); + wh.addValueToStorage(c(1), 'i_draw', v_bias); + wh.addValueToStorage(v(1), 'v_set', v_bias); + % Retrieve all stored output power measurements up to the current point - p_out = wh.getStoValue('p_out', wh.parameter.v_bias.values); - + p_opt_out = wh.getStoValue('p_out', wh.parameter.v_bias.values); + % Plot the measured output power in dBm versus the negative bias voltage - figure(90); - plot(-wh.parameter.v_bias.values(1:numel(p_out)), p_out, 'DisplayName', 'Measured Output Power in dBm'); - xlim([-max(params.v_bias), -min(params.v_bias)]); % Set x-axis limits + % Plot the output power converted from dBm to linear scale (Watts) + figure(91); + clf + hold on + plot(wh.parameter.v_bias.values(1:numel(p_opt_out)), db2pow(p_opt_out), 'DisplayName', 'Measured Output Power in Watts','LineWidth',2); + xlim([min(params.v_bias), max(params.v_bias)]); % Set x-axis limits grid on; % Enable grid for better readability + xlabel('applied voltage in V') + ylabel('measured optical power in dB') + end +save(['C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\HighSpeedMeasurement_2024\Messreihe_TFLN_MZM\',measurement_name],'wh'); + +p_opt_in = 0; %dbm % After completing the measurements, retrieve all output power data -p_out = wh.getStoValue('p_out', wh.parameter.v_bias.values); +p_opt_out = wh.getStoValue('p_out', wh.parameter.v_bias.values); +transmission = p_opt_out; +i_draw = wh.getStoValue('i_draw', wh.parameter.v_bias.values); +v_set = wh.getStoValue('v_set', wh.parameter.v_bias.values); + +power_draw = i_draw .* v_set .* 1e3; %watt +resistance = v_set ./ i_draw; +i_draw = i_draw.* 1e3; % Plot the output power converted from dBm to linear scale (Watts) +cols = linspecer(4); +c = 1; +linstyle = '-'; + figure(91); -plot(-wh.parameter.v_bias.values(1:numel(p_out)), db2pow(p_out), 'DisplayName', 'Measured Output Power in Watts'); -xlim([-max(params.v_bias), -min(params.v_bias)]); % Set x-axis limits +hold on +subplot(2,2,1) +hold on +plot(power_draw, transmission, 'DisplayName', measurement_name,'LineWidth',1,'Marker','none','Color',cols(c,:),'LineStyle',linstyle); +title('MZM Attenuation') +xlim([min(power_draw), max(power_draw)]); % Set x-axis limits grid on; % Enable grid for better readability +xlabel('DUT Power (mW)'); +ylabel('Transmission in dB'); +legend + +subplot(2,2,2) +hold on +plot(v_set, i_draw, 'DisplayName', measurement_name,'LineWidth',1,'Marker','none','Color',cols(c,:),'LineStyle',linstyle); +title('Voltage vs. Current') +xlim([min(v_set), max(v_set)]); % Set x-axis limits +grid on; % Enable grid for better readability +xlabel('Voltage (mA)'); +ylabel('Current (mA)'); +legend + +subplot(2,2,3) +hold on +plot(power_draw, db2pow(transmission), 'DisplayName', measurement_name,'LineWidth',1,'Marker','none','Color',cols(c,:),'LineStyle',linstyle); +title('Measured Output Power in Watt') +xlim([min(power_draw), max(power_draw)]); % Set x-axis limits +grid on; % Enable grid for better readability +xlabel('DUT Power (mW)'); +ylabel('Trabnsmission in %'); +legend + +subplot(2,2,4) +hold on +plot(power_draw, resistance, 'DisplayName', measurement_name,'LineWidth',1,'Marker','none','Color',cols(c,:),'LineStyle',linstyle); +title('Heater Resistance = U/I') +xlim([min(power_draw), max(power_draw)]); % Set x-axis limits +grid on; % Enable grid for better readability +xlabel('DUT Power (mW)'); +ylabel('Resistance in Ohm'); +ylim([455, 470]) +legend + +disp('Measurement Complete!') diff --git a/projects/Lab_2024/lab_sir_sweep.m b/projects/Lab_2024/lab_sir_sweep.m index 2bb3efc..ff1f0dd 100644 --- a/projects/Lab_2024/lab_sir_sweep.m +++ b/projects/Lab_2024/lab_sir_sweep.m @@ -1,13 +1,13 @@ -folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\sir_sweep_pam4\'; -experiment_name = 'PAM4_DB_encoded_10km_'; +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\sir_sweep_sd40\'; +experiment_name = 'PAM4_56_v2'; %%% SIR Sweep for MPI Experiment %%% params = struct; -params.vbias = [2.45]; -params.awg_vpp = [0.25]; -params.eq_mode = [4]; +params.vbias = [2.5]; +params.awg_vpp = [0.35]; +params.eq_mode = [2]; params.i_atten = [0:4:40]; wh = DataStorage(params); @@ -23,13 +23,13 @@ wh.addStorage("signals"); precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\standard_system\"; precomp_fn = "lab_mpi_setup_2"; precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active -precomp_amp_max = 5; +precomp_amp_max = 2; M = 4; pn_key = 2; usemrds = 0; -fsym = 92e9; -fdac = 92e9; +fsym = 68e9; +fdac = 256e9; awg_vpp = 0.35; fadc = 160e9; rrcalpha = 0.05; @@ -60,21 +60,6 @@ for eq_mode = wh.parameter.eq_mode.values sprintf('Progress: %d/%d\nEstimated time remaining: %.2f hours\nEstimated time remaining: %.2f hours', ... loopcnt, looptotal, estimatedTimeRemaining/60/60, estimatedTotalTime/60/60)); - try - - z - ber = wh2.getStoValue('ber',v_bias,awg_vpp,eq_mode,precomp_amp_max); - assert(~isempty(ber),'err') - rop = wh2.getStoValue('rop',v_bias,awg_vpp,eq_mode,precomp_amp_max); - assert(~isempty(rop)) - pd_in = wh2.getStoValue('pd_in',v_bias,awg_vpp,eq_mode,precomp_amp_max); - assert(~isempty(pd_in)) - M = wh2.getStoValue('m',v_bias,awg_vpp,eq_mode,precomp_amp_max); - assert(~isempty(M)) - - catch - - switch eq_mode @@ -93,6 +78,7 @@ for eq_mode = wh.parameter.eq_mode.values db_coding_approach = 0; db_precode = db_coding_approach || db_channel_approach; + case 3 ffe_only = 0; postfilter_approach = 0; @@ -122,8 +108,9 @@ for eq_mode = wh.parameter.eq_mode.values %%%%% Construct AWG and Scope Modules %%%%%% SCP = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc",'GSa_160',"channel",[1,0],"recordLen",2000000,"removeDC",1); - AWG = AwgKeysight("model","M8196A","fdac",fdac,"scaletodac",[1,1,1,1],"skews",[0,0,0,0],"voltages",[0,0,0,awg_vpp]); - A2S = Awg2Scope(AWG,SCP,[0,0,0,1]); + %AWG = AwgKeysight("model","M8196A","fdac",fdac,"scaletodac",[1,1,1,1],"skews",[0,0,0,0],"voltages",[0,0,0,awg_vpp]); + AWG = AwgKeysight("model","M8199B","fdac",fdac,"scaletodac",[1,1],"skews",[0,0],"voltages",[0,awg_vpp]); + A2S = Awg2Scope(AWG,SCP,[0,1,0,0]); %%%%% Symbol Generation %%%%%% Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); @@ -146,10 +133,12 @@ for eq_mode = wh.parameter.eq_mode.values %%%%% Resample to DAC rate %%%%%% Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + Digi_sig = Filter('filtdegree',5,"f_cutoff",1.1*(fsym/log2(M)),"fs",Digi_sig.fs,"filterType",filtertypes.gaussian).process(Digi_sig); + %%%%% Plot and Save Routine 1 %%%%%%%%%%%%%%%%%%%%%%%%% Digi_sig.spectrum("displayname","Normal Tx","fignum",10); - + if loopcnt == 1 save([folderpath,[experiment_name,'bits'],loop_name],"Bits"); save([folderpath,[experiment_name,'symbols'],loop_name],"Symbols"); @@ -157,7 +146,7 @@ for eq_mode = wh.parameter.eq_mode.values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% AWG --> Scope %%%%%% - [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); + [~,Scpe_sig,~,~] = A2S.process("signal2",Digi_sig); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); @@ -194,17 +183,43 @@ for eq_mode = wh.parameter.eq_mode.values % Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope'); voa.readvals(); - + pd_in = voa.power_state(2); s_pow = voa.power_state(3); i_pow = voa.power_state(4); - + sir = s_pow- i_pow; %%%%% EQUALIZE %%%%%% Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",50,"sps",2,"decide",0); - Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,7,7],"sps",2,"decide",1); + % Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,7,7],"sps",2,"decide",1); Eq = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + Eq = EQ("Ne",[50,0,0],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + + S = Scpe_sig.signal; + % recursion + N1 = 201; + % Initialize the moving sum for the first window + half_window = (N1 - 1) / 2; + moving_sum = sum(S(1:N1)); + + % Calculate the first element of R1 + cic_filtrd(1:half_window+1) = S(1:half_window+1) - (moving_sum / N1); + + % Loop over the signal and apply the recursive moving average subtraction + for n = (half_window+2):(length(S)-half_window) + % Update the moving sum by subtracting the oldest value and adding the new one + moving_sum = moving_sum - S(n-half_window-1) + S(n+half_window); + + % Calculate the new value of R1 + cic_filtrd(n) = S(n) - (moving_sum / N1); + end + cic_filtrd(n+1:length(S)) = S(n+1:length(S)) - (moving_sum / N1); + + Scpe_sig_ = Scpe_sig; + Scpe_sig_.signal = cic_filtrd'; + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -212,7 +227,6 @@ for eq_mode = wh.parameter.eq_mode.values EQ_sig.plot("fignum",50,"displayname",'After EQ'); - Rx_bits = PAMmapper(M,0).demap(EQ_sig); [~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); @@ -231,46 +245,62 @@ for eq_mode = wh.parameter.eq_mode.values elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% - [EQ_sig] = Eq.process(Scpe_sig,Symbols); + [EQ_sig, Noi] = Eq.process(Scpe_sig_,Symbols); - EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + % EQ_sig.plot("displayname",'After VNLE','fignum',90,'clear',1); - Noi = EQ_sig-Symbols; + % Quantization is too far from orig. symbols -> + % error psd is quite different + % Sym_ = PAMmapper(M,0).quantize(EQ_sig); + % Noi_ = Sym_-EQ_sig; + % Noi_.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:)); Rx_bits = PAMmapper(M,0).demap(EQ_sig); - [~,errors_bm,ber_ffe_only,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + [~,~,ber_vnle,~] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - nc = 2; - burg_coeff = arburg(Noi.signal,nc); + EQ_sig.normalize('mode','rms').spectrum('displayname','EQ Out','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc,:)); - EQ_sig = EQ_sig.filter(burg_coeff,1); + Noi.normalize('mode','rms').spectrum('displayname','Noise PSD optimal','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:)); - if 1 - Noi.spectrum('displayname','Noise PSD','fignum',123) - [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs); - h = h/max(abs(h)); - hold on - w_ = (w - Noi.fs/2); - plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); - end + for nc = 1:3 - if 0 - figure(53); - constellation = unique(Symbols.signal); - received = NaN(numel(constellation),length(Symbols)); - for lvl = 1:numel(constellation) - received(lvl,Symbols.signal==constellation(lvl)) = EQ_sig.signal(Symbols.signal==constellation(lvl)); + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig_filt = EQ_sig.filter(burg_coeff,1); + + % EQ_sig.spectrum("displayname","Signal Spectrum after Postfilter","fignum",1234); + + EQ_sig_mlse = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig_filt); + + % EQ_sig.spectrum("displayname","Signal Spectrum after MLSE","fignum",1234); + + if 0 + cols = linspecer(12); + + EQ_sig_filt.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+2,:)); + + % [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs); + % h = h/max(abs(h)); + % hold on + % w_ = (w - Noi.fs/2); + % figure(123) + % plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); + + [h,w] = freqz(1,burg_coeff,length(Noi),"whole"); + h = h/max(abs(h)); hold on - histogram(received(lvl,:),1000,"EdgeAlpha",0); + w_ = (w - pi); + plot(w_,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); + end + + Rx_bits = PAMmapper(M,0).demap(EQ_sig_mlse); + [~,errors_bm,ber(nc),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % disp(['BER: ',sprintf('%.1E',ber_mlse(i,j)),' - - ROP: ',num2str(patten(i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + end - EQ_sig = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); - - Rx_bits = PAMmapper(M,0).demap(EQ_sig); - [~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - disp(['FFE: ',sprintf('%.1E',ber_ffe_only),' -> PF -> MLSE: ',sprintf('%.1E',ber),' dB | PD_in: ',num2str(pd_in),' dBm']); + disp(['FFE: ',sprintf('%.1E',ber_vnle),' -> PF -> MLSE: ',sprintf('%.1E',ber(nc)),' dB | PD_in: ',num2str(pd_in),' dBm']); elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -301,8 +331,8 @@ for eq_mode = wh.parameter.eq_mode.values end - end + wh.addValueToStorage(ber_vnle,'vnle',v_bias,awg_vpp,eq_mode,i_atten); wh.addValueToStorage(ber,'ber',v_bias,awg_vpp,eq_mode,i_atten); wh.addValueToStorage(pd_in,'pd_in',v_bias,awg_vpp,eq_mode,i_atten); wh.addValueToStorage(Rx_bits,'signals',v_bias,awg_vpp,eq_mode,i_atten); @@ -312,13 +342,13 @@ for eq_mode = wh.parameter.eq_mode.values wh.addValueToStorage(M,'m',v_bias,awg_vpp,eq_mode,i_atten); showCurrentMeasurement('BER', ber, 'PD in', pd_in, 'PAM',M, 'Vbias', v_bias, 'AWG Vpp', awg_vpp, 'Precomp MaxAmp',precomp_amp_max); - + iterationTimes(loopcnt) = toc(iterationStartTime); averageTimePerIteration = mean(iterationTimes(1:loopcnt)); estimatedTotalTime = averageTimePerIteration * looptotal; estimatedTimeRemaining = estimatedTotalTime - sum(iterationTimes(1:loopcnt)); %autoArrangeFigures(3,3,2); - + end end end @@ -337,18 +367,12 @@ awg_vpp = wh.parameter.awg_vpp.values(1); eq_mode = wh.parameter.eq_mode.values(1); bers = wh.getStoValue('ber',v_bias,awg_vpp,eq_mode,i_atten_vals); - +sirs = wh.getStoValue('sir',v_bias,awg_vpp,eq_mode,i_atten_vals); figure(90); hold on; % Retain the plot so new points can be added without complete redraw % Plot the data and get the line handle -hLine = plot(i_atten_vals, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", experiment_name); - -% Customize the data tips -% Set labels for existing data tip rows -hLine.DataTipTemplate.DataTipRows(1).Label = 'Fsym'; -hLine.DataTipTemplate.DataTipRows(2).Label = 'BER'; -hLine.DataTipTemplate.DataTipRows(2).Format = '%.2e'; % Format BER as "3e-4" +plot(sirs, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", experiment_name); % Continue with the rest of your plot settings diff --git a/projects/standard_system/lab_mpi_setup_2.mat b/projects/standard_system/lab_mpi_setup_2.mat index d91e330..d399a78 100644 Binary files a/projects/standard_system/lab_mpi_setup_2.mat and b/projects/standard_system/lab_mpi_setup_2.mat differ diff --git a/wh.mat b/wh.mat deleted file mode 100644 index 320d18a..0000000 Binary files a/wh.mat and /dev/null differ