diff --git a/Classes/00_signals/Electricalsignal.m b/Classes/00_signals/Electricalsignal.m index 95beb14..7da0e69 100644 --- a/Classes/00_signals/Electricalsignal.m +++ b/Classes/00_signals/Electricalsignal.m @@ -3,7 +3,7 @@ classdef Electricalsignal < Signal % Detailed explanation goes here properties - fs + end methods @@ -26,7 +26,7 @@ classdef Electricalsignal < Signal end - function pow = power(obj) + function pow = power50(obj) % Power of an electrical signal R = 50; diff --git a/Classes/00_signals/Frame 1.png b/Classes/00_signals/Frame 1.png new file mode 100644 index 0000000..b1797d7 Binary files /dev/null and b/Classes/00_signals/Frame 1.png differ diff --git a/Classes/00_signals/Informationsignal.m b/Classes/00_signals/Informationsignal.m index 03ccc00..d098d2f 100644 --- a/Classes/00_signals/Informationsignal.m +++ b/Classes/00_signals/Informationsignal.m @@ -3,7 +3,7 @@ classdef Informationsignal < Signal % Detailed explanation goes here properties - fs + end methods @@ -26,11 +26,11 @@ classdef Informationsignal < Signal end - function pow = power(obj) - - pow = mean(abs(obj.signal.^2),"all") ; - - end + % function pow = power(obj) + % + % pow = mean(abs(obj.signal.^2),"all") ; + % + % end end diff --git a/Classes/00_signals/Opticalsignal.m b/Classes/00_signals/Opticalsignal.m index 7afc4b5..ee6137f 100644 --- a/Classes/00_signals/Opticalsignal.m +++ b/Classes/00_signals/Opticalsignal.m @@ -5,7 +5,6 @@ classdef Opticalsignal < Signal properties nase lambda - fs end @@ -52,14 +51,6 @@ classdef Opticalsignal < Signal end - - function pow = power(obj) - - pow = mean( abs(obj.signal).^2 ) ; - pow = pow2db(pow)+30; %dbm - - end - function cspr = cspr(obj) carrier_power_dbm = pow2db( abs(mean(obj.signal)).^2 )+30; % dB -> +30 -> dBm @@ -70,17 +61,15 @@ classdef Opticalsignal < Signal cspr_lin = pow2db(carr_lin/sign_lin); + c = abs(mean(obj.signal)).^2; + s = mean(abs(obj.signal).^2); + cspr_ = 10*log10(c / s); + cspr =carrier_power_dbm-signal_power_dbm; end - function papr = papr(obj) - average_power = obj.power; - peak_power = pow2db(max(abs(obj.signal.^2)))+30; - papr = peak_power - average_power; - - end end end diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 758bb96..bafbef7 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -5,14 +5,21 @@ classdef Signal properties signal logbook + fs end methods - function obj = Signal(signal) + function obj = Signal(signal,options) %SIGNAL Construct an instance of this class % Detailed explanation goes here + arguments + signal + options.fs = []; + end + obj.signal = signal; obj.signal = obj.signal; + obj.fs = options.fs; SignalType = []; TimeStamp = []; @@ -114,6 +121,64 @@ classdef Signal end + function plot(obj, options) + % signal to plot: obj.signal + % fsamp : obj.fs (e.g. 92e9 => 92 GHz) + % length: length(obj.signal) + + arguments + obj + options.fignum + options.displayname = []; + options.timeframe = 0; + end + + figure(options.fignum); % If figure does not exist, create new figure + + % 2) Plot into the figure handle found or created in one + t = (0:length(obj.signal)-1) / obj.fs; % time vector + + if options.timeframe ~= 0 + %only show a certain timeframe of signal + t = t(t I want to add the new - %spectum "onto" the existing plot to have a better comparison - if isempty(options.figurename) - fig = findall(groot, 'Type', 'figure', 'Name', 'power density'); - if isvalid(fig) - fig = get(fig); - ax = gca; - hold on - else - figure('name','power density'); - ax = gca; - end - else - fig = findall(groot, 'Type', 'figure', 'Name', options.figurename); - if isvalid(fig) - ax = fig.CurrentAxes; - hold on - else - figure('name',options.figurename); - ax = gca; - hold on - end - end + % spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname); + N = 2^(nextpow2(length(obj.signal))-6); + + + [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); + p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power" - %compute FFT of input - Fsignal = fft(obj.signal); - - %POWER spectral density (todo: toggle?) - psd = Fsignal.*conj(Fsignal); - - %Use only magnitude of FFT (which was complex) - psd = abs(psd); - - %Shift the spectrum to yield - psd = fftshift(psd); - - %divide by N - psd = psd/length(Fsignal); - - %smoothing - psd = smooth(psd,1000); - - psd_plot = 20*log10(psd); - - psd_plot(psd_plot<-120) = -120; - - - testParseval = 1; - if testParseval == 1 - E_FreqDomain = sum(psd); - %test parseval - E_TimeDomain = sum(abs(Fsignal.^2)); - - if isequal(round(E_FreqDomain,1),round(E_TimeDomain,1)) - %disp('Parseval is right!'); - else - disp('Parseval theorem is not right...'); - end - end - - if fsamp <= 1e+100 - %Frequency Axis - freq_vec = linspace(-fsamp/2,fsamp/2,length(psd)); - freq_vec = reshape(freq_vec,size(psd_plot)); - - - if ~isempty(options.displayname) - plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5,'DisplayName',options.displayname,'Parent',ax); - else - plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5,'Parent',ax); - end - - xlabel('Frequency [GHz]') - else - %Wavelength Axis - freq_vec = physconst('LightSpeed')*linspace(-fsamp/2,fsamp/2,length(psd))./((physconst('LightSpeed')/1550e-9)^2); - if ~isempty(options.displayname) - plot(freq_vec*1e9,psd_plot,'Linewidth',0.5,'DisplayName',options.displayname,'Parent',ax); - else - plot(freq_vec*1e9,psd_plot,'Linewidth',0.5,'Parent',ax); - end - - xlabel('Wavelength [nm]') - end - - - ylabel('Magnitude [dB]') + figure(options.fignum); % If figure does not exist, create new figure + hold on + plot(w.*1e-9,p_dbm,'DisplayName',options.displayname,'LineWidth',1); + xlabel("Frequency in GHz"); + %ylabel("Power/frequency (dB/Hz)"); + ylabel("Power (dBm)"); + xlim([-obj.fs/2 obj.fs/2].*1e-9) + edgetick = 2^(nextpow2(obj.fs*1e-9)); + xticks([-edgetick:16:edgetick]); + xlim([-244, 244]) + ylim([-120,-0]); + yticks([-200:10:10]); legend - grid minor; end + %% Power of signal + function pow = power(obj,options) + + arguments + obj + options.unit power_notation = power_notation.dBm + end + + pow = mean(abs(obj.signal).^2); + + switch options.unit + case power_notation.dBm + if isa(obj,'Electricalsignal') + pow = pow / 50; + end + pow = 10*log10(pow)+30; %dbm + + case power_notation.mW + pow = pow .* 1e3; %mW + + case power_notation.W + %pow = pow % Watt + + end + + end + + %% Peak Power of Signal + function pow_pk = power_peak(obj,options) + + arguments + obj + options.unit power_notation = power_notation.dBm + end + + pow_pk = max(abs(obj.signal).^2); % dBm + + switch options.unit + case power_notation.dBm + pow_pk = pow2db(pow_pk)+30; %dbm + case power_notation.mW + pow_pk = pow_pk .* 1e3; %mW + case power_notation.W + %pow = pow % Watt + end + + end + + %% PAPR of signal + function papr = papr_lin(obj) + %PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power) + % divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor. + % papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2; + + papr = obj.power_peak("unit",power_notation.W) / obj.power("unit",power_notation.W); %linear + + end + + %% PAPR of signal + function papr_db = papr_db(obj) + %PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power) + % divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor. + % papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2; + + papr = obj.power_peak("unit",power_notation.W) / obj.power("unit",power_notation.W); + papr_db = 10*log10(papr); + + average_power = obj.power; + peak_power = obj.power_peak; + papr_db = peak_power - average_power; %db + + end %% function obj = normalize(obj,options) @@ -366,33 +434,310 @@ classdef Signal end - function eye(obj,fsym) + function er = extinctionratio(obj,fsym,M) + histpoints = 1024; %% verticale resolution + histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven + histpoints_horizontal = 512; %% horizontal resolution + hist_data=zeros(histpoints,histpoints_horizontal ); %% initilize eye diagram - disp('h'); - - fsig = obj.fs; - q = fsig/fsym; - - if q > 10 && isinteger(q) - sig = (obj.signal); + if isa(obj,'Opticalsignal') + sig = abs(obj.signal).^2; + elseif isa(obj,'Electricalsignal') + sig = obj.signal; else - sig = (obj.resample("fs_in",fsig,"fs_out",fsym*30).signal); - q = 10; + sig = obj.signal; end - figure() + x = (sig); %% make input signal rea)l + + x = resample(x,fsym*histpoints_horizontal/2,obj.fs); %% up sample to original fsym rate + + if mod(length(x),2)==1 %% if the signal lenght is not divisible by 2 (symbols displayed in the eye diagram are 2) remove last symbol + x = x(1:end-1); + end + + eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols) + + maxA = max(sig(100:end-100)); + minA = min(sig(100:end-100)); + + difference= maxA-minA; + + data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1; + + for n=1:size(data_ind_y,1) + nn=histcounts(data_ind_y(n,:),1:histpoints+1); + hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-( + end + + plot_data = 20*log10(hist_data); + plot_data(plot_data==-Inf) = 0; + + maxall = 0; + + for l = 1:size(plot_data,2) + [maxpk_,pos_] = max(plot_data(:,l)); + if maxpk_ > maxall + maxall = maxpk_; + posxall = l; + posyall = pos_; + end + end + + hist_interest = plot_data(:,posxall); + hist_interest_smoth = smooth(hist_interest,20); + [pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30); + + for i = 1:numel(loc) + ppeak(i) = maxA - (difference/histpoints*loc(i)); + end + + + + if isa(obj,'Opticalsignal') + er=10*log10(ppeak(1)/ppeak(end)); + elseif isa(obj,'Electricalsignal') + if mean([ppeak(1),ppeak(end)]) < 1e-2 + disp("No Extiction Ration for Bipolar Electrical Signal. Calculating Outer OMA instead...") + er=max(ppeak)-min(ppeak); + else + er=10*log10(ppeak(1)/ppeak(end)); + end + else + er=10*log10(ppeak(1)/ppeak(end)); + end + + if 0 + findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30); + end + + + end + + function eye(obj,fsym,M) + + mode = 1; + + histpoints = 1024; %% verticale resolution + histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven + histpoints_horizontal = 512; %% horizontal resolution + hist_data=zeros(histpoints,histpoints_horizontal ); %% initilize eye diagram + + if isa(obj,'Opticalsignal') + sig = abs(obj.signal).^2; + elseif isa(obj,'Electricalsignal') + sig = obj.signal; + else + sig = obj.signal; + end + + x = (sig); %% make input signal rea)l + + x = resample(x,fsym*histpoints_horizontal/2,obj.fs); %% up sample to original fsym rate + + if mod(length(x),2)==1 %% if the signal lenght is not divisible by 2 (symbols displayed in the eye diagram are 2) remove last symbol + x = x(1:end-1); + end + + eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols) + + figure(922) clf - cursor = 200*q; + if mode == 2 + % generate "intuitive eye diagram" by drawing lines on top over + % each other; only draw 1000 lines, otherwise the plot is too + % crowded + + + col = cbrewer2('Set1',2); + for n=1:1000 + hold on + plot(eye_mat(:,n),'LineStyle',':','LineWidth',0.1,'Color',col(2,:)); + end + ylabel('Amplitude of Signal'); + + elseif mode == 1 + % generate eye diagram using histogram + + maxA = max(sig(100:end-100)); + minA = min(sig(100:end-100)); + + + + difference= maxA-minA; + + data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1; + + for n=1:size(data_ind_y,1) + nn=histcounts(data_ind_y(n,:),1:histpoints+1); + hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-( + end + + plot_data = 20*log10(hist_data); + plot_data(plot_data==-Inf) = 0; + + imagesc(plot_data); + + + + % beautify + colormap(cbrewer2("Blues",4096)); + + + if isa(obj,'Opticalsignal') + title("Optical Eye") + ylabel("Power in mW"); + y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,16)); + min_ = min(abs(obj.signal(100:end-100)).^2); + max_ = abs(max(obj.signal(100:end-100)).^2); + elseif isa(obj,'Electricalsignal') + title("Electrical Eye") + ylabel("Voltage in V"); + y_tickstring = string(linspace(maxA,minA,16)); + min_ = min(obj.signal(100:end-100)); + max_ = abs(max(obj.signal(100:end-100))); + else + title("Digital Eye") + ylabel("Digital Signal Amplitude"); + y_tickstring = string(linspace(maxA,minA,16)); + min_ = min(obj.signal(100:end-100)); + max_ = abs(max(obj.signal(100:end-100))); + end + + % add information + + if 1 + + pwr_dbm = round(obj.power,3); + pwr_lin = obj.power("unit",power_notation.W); + + papr_ = obj.papr_lin;%round(papr(obj.signal.^2),3); + + yline( histpoints-(pwr_lin - minA)/difference*histpoints ); + yline( histpoints-(min_ - minA)/difference*histpoints ); + yline( histpoints-(max_ - minA)/difference*histpoints ); + + maxall = 0; + for l = 1:size(plot_data,2) + [maxpk_,pos_] = max(plot_data(:,l)); + if maxpk_ > maxall + maxall = maxpk_; + posxall = l; + posyall = pos_; + end + end + + hold on + xline(posxall) + + hist_interest = plot_data(:,posxall); + hist_interest_smoth = smooth(hist_interest,20); + a = scatter(hist_interest_smoth+posxall,1:length(hist_interest_smoth),4,'.','MarkerEdgeColor','red'); + + [pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30); + scatter(posxall,loc,'red','Marker','x','LineWidth',2); + yline(loc,'Color','red','LineWidth',1,'LineStyle',':'); + + for i = 1:numel(loc) + ppeak(i) = maxA - (difference/histpoints*loc(i)); + end + + oma = false; + if isa(obj,'Opticalsignal') + er=10*log10(ppeak(1)/ppeak(end)); + elseif isa(obj,'Electricalsignal') + if mean([ppeak(1),ppeak(end)]) < 1e-2 + oma = true; + er=max(ppeak)-min(ppeak); + else + er=10*log10(ppeak(1)/ppeak(end)); + end + else + er=10*log10(ppeak(1)/ppeak(end)); + end + + % Define properties + boxPosition = [0.15 0.86 0.2 0.05]; % Position for the first box [x y width height] + boxColor = [0.9 0.9 0.9]; % Light grey background color + boxEdgeColor = 'k'; % Black edge color + boxLineStyle = '--'; % Dashed line style + boxFontWeight = 'bold'; % Bold font + + % Create first annotation box for Power + annotation('textbox', boxPosition, ... + 'String', ['Power: ',num2str(pwr_dbm),' dBm'], ... + 'BackgroundColor', boxColor, ... + 'EdgeColor', boxEdgeColor, ... + 'LineStyle', boxLineStyle, ... + 'FontWeight', boxFontWeight, ... + 'HorizontalAlignment', 'center'); + + % Adjust position for the second box (slightly to the right) + boxPosition = [0.37 0.86 0.2 0.05]; % Adjusted position + + % Create second annotation box for PAPR + annotation('textbox', boxPosition, ... + 'String', ['PAPR(lin):',num2str(papr_),''], ... + 'BackgroundColor', boxColor, ... + 'EdgeColor', boxEdgeColor, ... + 'LineStyle', boxLineStyle, ... + 'FontWeight', boxFontWeight, ... + 'HorizontalAlignment', 'center'); + + % Adjust position for the third box (slightly to the right) + boxPosition = [0.59 0.86 0.2 0.05]; % Adjusted position + + % Create third annotation box for Vmax + if ~oma + thirdboxstring = ['ER (db):',num2str(er),' dB']; + else + thirdboxstring = ['OMA outer:',num2str(er),' V']; + end + + annotation('textbox', boxPosition, ... + 'String',thirdboxstring , ... + 'BackgroundColor', boxColor, ... + 'EdgeColor', boxEdgeColor, ... + 'LineStyle', boxLineStyle, ... + 'FontWeight', boxFontWeight, ... + 'HorizontalAlignment', 'center'); + + + yticks(linspace(0,histpoints,16)); + yticklabels(y_tickstring); + + end + + + - for s = 1:700 - plot(sig(cursor-q:cursor+q),'Color','black','LineWidth',0.1,'LineStyle','-'); - hold on - cursor=cursor+q; - s=s+1; end - ylim([-3 3]); - + % disp('h'); + % + % fsig = obj.fs; + % q = fsig/fsym; + % + % if q > 10 && isinteger(q) + % sig = (obj.signal); + % else + % sig = (obj.resample("fs_in",fsig,"fs_out",fsym*30).signal); + % q = 10; + % end + % + % figure() + % clf + % cursor = 200*q; + % + % for s = 1:700 + % plot(sig(cursor-q:cursor+q),'Color','black','LineWidth',0.1,'LineStyle','-'); + % hold on + % cursor=cursor+q; + % s=s+1; + % end + % + % ylim([-3 3]); + end diff --git a/Classes/01_transmit/AWG.m b/Classes/01_transmit/AWG.m index c82b5b3..79c5469 100644 --- a/Classes/01_transmit/AWG.m +++ b/Classes/01_transmit/AWG.m @@ -46,7 +46,7 @@ classdef AWG < handle options.lpf_type = 0; options.f_cutoff = 32e9; options.H_lpf Filter - + end fn = fieldnames(options); @@ -60,26 +60,33 @@ classdef AWG < handle function signalclass_out = process(obj,signalclass_in) + + % 0) START AWG SIMULATION + len_in = length(signalclass_in.signal); + % disp(['Power Raw Input: ', num2str(mean(abs(real(signalclass_in.signal).^2)))]); + + % 1) RESAMPLE IF NESSECARY 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); + % min_ = min(signalclass_in.signal); + % max_ = max(signalclass_in.signal); + signalclass_in = signalclass_in.resample("fs_in",signalclass_in.fs,"fs_out",obj.fdac,"n",10,"beta",5); + % signalclass_in.signal = clip(signalclass_in.signal,-1.3414,1.3414); + % signalclass_in.signal = softclip(signalclass_in.signal,7,1.3414); + % signalclass_in.plot("displayname",'resampled and clipped','fignum',1); end - - % 1-3. actual processing of the signal (normalize->quantize->sample hold) + + % disp(['Power after resamp: ', num2str(mean(abs(real(signalclass_in.signal).^2)))]); + + % 4) PROCESSING of the signal (normalize->quantize->sample hold) 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(13,"dBm"); - - % 4. Apply LPF on the signal if obj.lpf_active if isa(obj.H_lpf,'Filter') @@ -92,9 +99,9 @@ classdef AWG < handle signalclass_in = lpf.process(signalclass_in); end end - - disp(['AWG output power: ',num2str(signalclass_in.power),' dBm']); - + + % disp(['AWG output power: ',num2str(signalclass_in.power),' dBm']); + % append to logbook current_class = class(obj); lbdesc = ['AWG ', current_class , '// k_over:',num2str(obj.kover),'. f_dac:',num2str(obj.fdac*1e-9),'GHz. Resolution:',num2str(obj.bit_resolution),' bits.']; @@ -125,9 +132,57 @@ classdef AWG < handle obj.signal_length = length(data_in); + %%%%%%%%% PRECOMP SINC ROLLOFF %%%%%%%%% + if 1 + % X: design FIR filter for sinc precomp + ntaps = 13; + npts = 32; + % least-squares FIR design + fmax = obj.fdac*0.5; + ff = linspace(0,fmax,npts); + hsinc = sin(pi*ff/obj.fdac)./(pi*ff/obj.fdac + eps); % transfer function of sample and hold DAC + hsinc(1) = 1; + h_goal= 1./hsinc; % goal function + f = 2.*ff./obj.fdac; %vector between 0 and 1, where 1 is nyquist is fsamp/2 + b = firls(ntaps-1,f,h_goal); + data_in = conv(data_in,b,"same"); + end + + if 0 + %compare different fir construction methods in matlab + b_fir = fir2(ntaps-1,f,h_goal); + b_pm = firpm(ntaps-1,f,h_goal); + b = firls(ntaps-1,f,h_goal); + + figure(111) + freqz(b,1,[],obj.fdac); + hold on + freqz(b_fir,1,[],obj.fdac); + hold on + freqz(b_pm,1,[],obj.fdac); + plot(f.*obj.fdac./2.*1e-9,20*log10(h_goal),'Marker','o'); + legend('firls','fir2','firpm','target') + end + + if 0 + % design filter as inverse of the goal transfer function + freq_vec = linspace(-obj.fdac/2,obj.fdac/2,length(data_in)); + h = sin(pi*freq_vec/obj.fdac)./(pi*freq_vec/obj.fdac + eps); + max_amp_db = 45; + max_amp_lin = 10^(-max_amp_db/20); + p = find(h0 @@ -144,9 +202,9 @@ classdef AWG < handle else elec_out = data_in; end - - % 2. Sample and hold + repeat (data_out: 1xsignal length) + %%%%%%%%% Sample and hold %%%%%%%%% + % X. Sample and hold + repeat (data_out: 1xsignal length) if obj.upsampling_method == 1 % just use matlab function elec_out = resample(elec_out,obj.kover,1); @@ -157,8 +215,8 @@ classdef AWG < handle else error('chosen upsampling method not implemented?'); end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 3. Add skew (not implemented so far) if obj.skew_active elec_out = obj.skew(elec_out); diff --git a/Classes/01_transmit/M8196A.m b/Classes/01_transmit/M8196A.m index 8a0d7f8..abe820c 100644 --- a/Classes/01_transmit/M8196A.m +++ b/Classes/01_transmit/M8196A.m @@ -5,26 +5,22 @@ classdef M8196A < AWG methods - function obj = M8196A() + function obj = M8196A(options) %This is a ready to use AWG simulation that resembles the %properties of the Keysighe M8196A % M8196A (92GBd) https://www.keysight.com/us/en/product/M8196A/92-gsa-s-arbitrary-waveform-generators.html arguments - %options.bla = 1; + options.kover = 8; end - obj = obj@AWG(); - obj.dac_max = 0.5; - obj.dac_min = -.5; - obj.bit_resolution = 5.5; + dac_max = 0.4; + dac_min = -0.4; - obj.fdac = 92e9; - - obj.lpf_active = 1; - - obj.f_cutoff = 32e9; - obj.lpf_type = filtertypes.gaussian; - obj.H_lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type); + fdac = 92e9; + + Lp_awg = Filter('filtdegree',4,"f_cutoff",32e9,"fs",fdac*options.kover,"filterType",filtertypes.butterworth,"active",true); + obj = obj@AWG("fdac",fdac,"dac_min",dac_min,"dac_max",dac_max,"lpf_active",1,"H_lpf",Lp_awg,"kover",options.kover,... + "bit_resolution",5.5,"normalize2dac",1,"upsampling_method","samplehold"); end diff --git a/Classes/01_transmit/M8199A.m b/Classes/01_transmit/M8199A.m new file mode 100644 index 0000000..e46121c --- /dev/null +++ b/Classes/01_transmit/M8199A.m @@ -0,0 +1,32 @@ +classdef M8199A < AWG + properties + + end + + methods + + function obj = M8199A(options) + %This is a ready to use AWG simulation that resembles the + %properties of the Keysight M8199A + + arguments + options.kover = 8; + end + + dac_max = 0.3; + dac_min = -0.3; + + fdac = 256e9; + + Lp_awg = Filter('filtdegree',4,"f_cutoff",65e9,"fs",fdac*options.kover,"filterType",filtertypes.butterworth,"active",true); + + obj = obj@AWG("fdac",fdac,"dac_min",dac_min,"dac_max",dac_max,"lpf_active",1,"H_lpf",Lp_awg,"kover",options.kover,... + "bit_resolution",5.5,"normalize2dac",1,"upsampling_method","samplehold"); + + end + + + end + + +end \ No newline at end of file diff --git a/Classes/01_transmit/M8199B.m b/Classes/01_transmit/M8199B.m index 448e692..2519453 100644 --- a/Classes/01_transmit/M8199B.m +++ b/Classes/01_transmit/M8199B.m @@ -5,25 +5,22 @@ classdef M8199B < AWG methods - function obj = M8199B() + function obj = M8199B(options) %This is a ready to use AWG simulation that resembles the %properties of the Keysighe M8199B https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf arguments - %options.bla = 1; + options.kover = 8; end - obj = obj@AWG(); - obj.dac_max = 0.5; - obj.dac_min = -.5; - obj.bit_resolution = 5.5; + dac_max = 0.6; + dac_min = -0.6; - obj.fdac = 256e9; + fdac = 256e9; + + Lp_awg = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*options.kover,"filterType",filtertypes.butterworth,"active",true); - obj.lpf_active = 1; - - obj.f_cutoff = 80e9; - obj.lpf_type = filtertypes.gaussian; - obj.H_lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type); + obj = obj@AWG("fdac",fdac,"dac_min",dac_min,"dac_max",dac_max,"lpf_active",1,"H_lpf",Lp_awg,"kover",options.kover,... + "bit_resolution",5.5,"normalize2dac",1,"upsampling_method","samplehold"); end diff --git a/Classes/01_transmit/PAMsource.m b/Classes/01_transmit/PAMsource.m new file mode 100644 index 0000000..ec7c55e --- /dev/null +++ b/Classes/01_transmit/PAMsource.m @@ -0,0 +1,117 @@ +classdef PAMsource + %NAME Summary of this class goes here + % Detailed explanation goes here + + properties(Access=public) + order + useprbs + M + fsym + randkey + + applypulseform + pulseformer + + fs_out + + applyclipping + clipfactor + end + + methods (Access=public) + function obj = PAMsource(options) + %NAME Construct an instance of this class + % Detailed explanation goes here + + arguments + options.order = 16; + options.useprbs = true; + options.M = true + options.fsym = 112e9; + options.randkey = 0; + + options.applypulseform = 1; + options.pulseformer Pulseformer; + + options.fs_out ; + + options.applyclipping = 0; + options.clipfactor = 10; + end + + + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + + if boolean(obj.applypulseform) && isempty(obj.pulseformer) + warning('No Pulseformer given. Proceeding with RRC and alpha 0.05'); + options.pulseformer = Pulseformer("fsym",obj.fsym,"fdac",obj.fs_out,"pulse","rrc","pulselength",16,"rrcalpha",0.05); + end + + + + end + + function [digi_sig,symbols,bits] = process(obj) + + %%%%% PRBS Generation in correct shape for Modulation Format %%%%%% + O = obj.order; %order of prbs + N = 2^(O-1); %length of prbs + [~,seed] = prbs(O,1); %initialize first seed of prbs + bitpattern=[]; + + if obj.useprbs + for i = 1:log2(obj.M) + [bitpattern(:,i),seed] = prbs(O,N,seed); + end + else + s = RandStream('twister','Seed',obj.randkey); + for i = 1:log2(obj.M) + bitpattern(:,i) = randi(s,[0 1], N, 1); + end + end + + if obj.M == 6 + bitpattern = reshape(bitpattern,[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); + end + + bits = Informationsignal(bitpattern); + + symbols = PAMmapper(obj.M,0).map(bits); + symbols.fs = obj.fsym; + + if obj.applyclipping + sym_min = min(symbols.signal); + sym_max = max(symbols.signal); + end + + %%%%% Pulseforming %%%%%% + if obj.applypulseform + digi_sig = obj.pulseformer.process(symbols); + else + digi_sig = symbols; + end + + %%%%% Resample to f DAC %%%%%% + digi_sig = digi_sig.resample("fs_in",digi_sig.fs,"fs_out",obj.fs_out,"n",10,"beta",5); + + %%%%% Hard clip digital signal to PAM range before DAC %%%%%% + if obj.applyclipping + digi_sig.signal = clip(digi_sig.signal , sym_min * obj.clipfactor , sym_max * obj.clipfactor); + end + + % append to logbook + lbdesc = ['Generated PAM Signal']; + digi_sig = digi_sig.logbookentry(lbdesc); + + end + + + end + +end diff --git a/Classes/01_transmit/Pulseformer.m b/Classes/01_transmit/Pulseformer.m index a3fa2cd..a306ef8 100644 --- a/Classes/01_transmit/Pulseformer.m +++ b/Classes/01_transmit/Pulseformer.m @@ -123,9 +123,9 @@ classdef Pulseformer %data_out_ = upfirdn(data_in,h,up,dn); % % %cut signal, which is longer due to fir filter - st = round(up/dn*racos_len/2); %we need to cut y_out -% en = round(st + (length(data_in)*up/dn) -1); -% data_out = data_out(st:en); + % st = round(up/dn*racos_len/2); %we need to cut y_out + % en = round(st + (length(data_in)*up/dn) -1); + % data_out = data_out(st:en); %scaling?! see pulsef module line 696 % scale = max(max([abs(real(data_out)) abs(imag(data_out))])); %find max value from real and imag part @@ -134,8 +134,8 @@ classdef Pulseformer data_out = data_out'; %Check output integrity - if round(up/dn * length(data_in)) ~= length(data_out) - %warning('Check signal length after pulse shaping'); + if abs(round(up/dn * length(data_in)) - length(data_out)) > 4 + warning('Check signal length after pulse shaping'); %disp('Check signal length after pulse shaping'); end diff --git a/Classes/01_transmit/Signalgenerator.m b/Classes/01_transmit/Signalgenerator.m new file mode 100644 index 0000000..6de1214 --- /dev/null +++ b/Classes/01_transmit/Signalgenerator.m @@ -0,0 +1,87 @@ +classdef Signalgenerator + %NAME Summary of this class goes here + % Detailed explanation goes here + + properties(Access=public) + form + length + fs + fsig + + end + + methods (Access=public) + function obj = Signalgenerator(options) + %NAME Construct an instance of this class + % Detailed explanation goes here + + arguments + options.form signalform = signalform.sine + options.length double = 1024 + options.fs double = 1000 %Hz sampling + options.fsig double = 50 % Hz fundamental frex e.g. of the sine or sawtooth + end + + % + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + + end + + function signalclass_out = process(obj) + + % actual processing of the signal (steps 1. - 3.) + signal = obj.build_signal(); + signalclass_out = Informationsignal(signal,"fs",obj.fs); + + % append to logbook + lbdesc = ['Signalgenerator: ',char(obj.form), '']; + signalclass_out = signalclass_out.logbookentry(lbdesc); + + end + + function signal = build_signal(obj) + %METHOD1 Summary of this method goes here + % Detailed explanation goes here + arguments(Input) + obj + end + + arguments(Output) + signal double + end + + switch obj.form + case signalform.sine + % Parameters + T = 1/obj.fs; % Sampling period (seconds per sample) + L = obj.length; % Length of signal (number of samples) + t = (0:L-1)*T; % Time vector + + % Sine wave parameters + f = obj.fsig; % Frequency of the sine wave (Hz) + A = 1; % Amplitude of the sine wave + + % Generate sine wave + signal = A * sin(2*pi*f*t); + case signalform.noise + + end + + + + end + + end + + methods (Access=private) + % Cant be seen from outside! So put all your functions here that can/ + % shall not be called from outside + + + end +end diff --git a/Classes/02_etc/Filter.m b/Classes/02_etc/Filter.m index 05ae2fb..6145e5b 100644 --- a/Classes/02_etc/Filter.m +++ b/Classes/02_etc/Filter.m @@ -3,6 +3,8 @@ classdef Filter < handle % Detailed explanation goes here properties + active + H filterType f_cutoff @@ -22,6 +24,7 @@ classdef Filter < handle %FILTER Construct an instance of this class % Detailed explanation goes here arguments + options.active logical = true; options.filterType filtertypes = filtertypes.bessel_inp ; options.f_cutoff = 0; options.fs = 0; @@ -47,7 +50,10 @@ classdef Filter < handle end function signal_out = process(obj,signal_in) - + if ~obj.active + signal_out = signal_in; + return + end if isa(signal_in,'Signal') % actual processing of the signal @@ -284,10 +290,12 @@ classdef Filter < handle title('Magnitude') hold on p = plot(frex,filt,'LineWidth',1,'LineStyle','-','DisplayName',[cur_module_name,': ',filter_desc, '; 3/6/10 dB: ', num2str(threedB),'/ ',num2str(sixdB),'/ ',num2str(ninedB)]); - + + if 0 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); + end xline(fcut_,'LineStyle',':','LineWidth',1,'HandleVisibility','off','Color',p.Color); yline([-3, -6, -9],'LineStyle',':','LineWidth',1,'HandleVisibility','off'); diff --git a/Classes/02_optical/Fiber.m b/Classes/02_optical/Fiber.m index eb0c67c..55f8ec2 100644 --- a/Classes/02_optical/Fiber.m +++ b/Classes/02_optical/Fiber.m @@ -49,7 +49,7 @@ classdef Fiber function signalclass_out = process(obj,signalclass_in) % actual processing of the signal (steps 1. - 3.) - signalclass_in.signal = obj.process_(signalclass_in.signal); + signalclass_in = obj.process_(signalclass_in); % append to logbook lbdesc = 'Fiber '; @@ -60,27 +60,45 @@ classdef Fiber end - function opt_out = process_(obj,opt_in) + function opt_sig = process_(obj,opt_sig) %METHOD1 Summary of this method goes here % Detailed explanation goes here - obj.b2 = -obj.D*obj.lambda0^2/(2*pi*Constant.LightSpeed); - obj.b3 = ((obj.lambda0.^2/(2*pi*Constant.LightSpeed)).^2*obj.Dslope); - obj.alpha_lin = obj.alpha/10*log(10)/1000; + signal = opt_sig.signal; + lambda_signal = opt_sig.lambda; - N = length(opt_in); + obj.D = obj.D + (lambda_signal-obj.lambda0)*obj.Dslope; + + obj.b2 = -obj.D*lambda_signal^2/(2*pi*Constant.LightSpeed); + obj.b3 = ((lambda_signal.^2/(2*pi*Constant.LightSpeed)).^2*obj.Dslope); + obj.alpha_lin = obj.alpha/10*log(10)/1000; + + N = length(signal); faxis = linspace(-obj.fsimu/2,obj.fsimu/2,N+1); faxis = ifftshift(faxis(:,1:end-1)); faxis = faxis'; - + obj.linstep = -obj.alpha_lin/2 - 2*1j*pi^2*obj.b2*faxis.^2 - 4/3*1j*pi^3*obj.b3*faxis.^3; - if obj.gamma ~= 0 - opt_out = obj.NLSE(opt_in); - else - opt_out = ifft(fft(opt_in).*exp(obj.linstep*obj.fiber_length)); % only one linear step + if 0 + H = exp((obj.linstep)*obj.fiber_length); + + figure(222) + hold on + plot(faxis.*1e-9,abs(real(Y)),'LineStyle','-','DisplayName','Abs(real) part of complex TF'); + xlabel("Frequency in GHz") + ylabel("$ R|(H(\omega, L))|$") end - %attenuate nase + + if obj.gamma ~= 0 + opt_out = obj.NLSE(signal); + else + opt_out = ifft( fft(signal) .* exp(obj.linstep*obj.fiber_length) ); % only one linear step + end + + %TODO: attenuate nase ... + + opt_sig.signal = opt_out; end diff --git a/Classes/03_receive/Photodiode.m b/Classes/03_receive/Photodiode.m index 9b02436..ecf9cac 100644 --- a/Classes/03_receive/Photodiode.m +++ b/Classes/03_receive/Photodiode.m @@ -41,6 +41,8 @@ classdef Photodiode % cast the inform. signal to electrical signal [signalclass_in, nase, lambda] = Electricalsignal(signalclass_in,"fs",obj.fsimu,"logbook",signalclass_in.logbook); + + % append to logbook lbdesc = ['Photo Diode ']; signalclass_in = signalclass_in.logbookentry(lbdesc); @@ -71,8 +73,8 @@ classdef Photodiode % Thermal Noise therm_current_psd = (2 * k * T / R ) ; %squared - 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 + Bw = obj.fsimu; + therm_noise_pow = therm_current_psd * Bw; %squared therm_noise = sqrt(therm_noise_pow) .* randn(obj.randomstream,size(yout,1),1); diff --git a/Classes/04_DSP/A1_scheme.m b/Classes/04_DSP/A1_scheme.m new file mode 100644 index 0000000..0a76c1c --- /dev/null +++ b/Classes/04_DSP/A1_scheme.m @@ -0,0 +1,23 @@ +classdef A1_scheme + %A1_SCHEME Summary of this class goes here + % Detailed explanation goes here + + properties + Property1 + end + + methods + function obj = A1_scheme(inputArg1,inputArg2) + %A1_SCHEME Construct an instance of this class + % Detailed explanation goes here + obj.Property1 = inputArg1 + inputArg2; + end + + function outputArg = method1(obj,inputArg) + %METHOD1 Summary of this method goes here + % Detailed explanation goes here + outputArg = obj.Property1 + inputArg; + end + end +end + diff --git a/Classes/04_DSP/EQ_silas.m b/Classes/04_DSP/EQ_silas.m index d20f438..42470af 100644 --- a/Classes/04_DSP/EQ_silas.m +++ b/Classes/04_DSP/EQ_silas.m @@ -170,23 +170,23 @@ classdef EQ_silas < handle for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength m = m+1; dc_cnt = dc_cnt+1; - + %get Sigal input vectors with correct length for VNLE x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).'; - + x_in_vnle_format = obj.calcVNLENonlinVecs(x_in_block,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); - + %get Reference input vectors with correct length for VNLE d_block = obj.d(obj.Nb(1)-obj.delay+m-2:-1:m-obj.delay-1).'; d_vnle_format = obj.calcVNLENonlinVecs(d_block,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); - + obj.e_ffe = obj.e.' * x_in_vnle_format; - + obj.e_dfe = obj.b.' * d_vnle_format; - + % Calculate the Error obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); - + if obj.mu_ffe_train ~= 0 %update FFE coefficients with LMS obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train; @@ -194,18 +194,17 @@ classdef EQ_silas < handle %update FFE coefficients with NLMS obj.e = obj.e - obj.error*x_in_vnle_format/(x_in_vnle_format.'*x_in_vnle_format); end - + %update DFE coefficients with LMS obj.b = obj.b + obj.mu_dfe_train*obj.error*d_vnle_format; - + %update DC error dc_block(dc_cnt) = obj.error .* obj.mu_dc_train; - + if dc_cnt == obj.eq_parallelization_blocklength obj.e_dc = obj.e_dc - mean(dc_block(dc_cnt)); dc_cnt = 0; end - end @@ -227,42 +226,29 @@ classdef EQ_silas < handle dc_cnt = 0; mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe - ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe - ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)... %3rd order ffe - ones(1,sum(obj.Cb))*obj.mu_dfe_dd]); %all order dfe - - mu_ffe = [ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe - ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe - ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)]; - - mu_dfe = ones(1,sum(obj.Cb))*obj.mu_dfe_dd; + ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe + ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)... %3rd order ffe + ones(1,sum(obj.Cb))*obj.mu_dfe_dd]); %all order dfe y = zeros(1,floor(obj.x_length/obj.sps)); d_feedback = zeros(obj.Cb(1),1); d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); - d_hat = zeros(obj.x_length,1); - lvl_err = NaN(obj.x_length,numel(obj.d_constellation)); - lvl_err_mov = NaN(100,numel(obj.d_constellation)); + d_hat = NaN(length(obj.d),numel(obj.d_constellation)); + lvl_err_1 = NaN(length(obj.d),numel(obj.d_constellation)); + lvl_err_2 = NaN(length(obj.d),numel(obj.d_constellation)); + subtracted_error =NaN(length(obj.d),numel(obj.d_constellation)); + y_1= NaN(length(obj.d),numel(obj.d_constellation)); + y_2= NaN(length(obj.d),numel(obj.d_constellation)); + lvl_err_mov = NaN(obj.eq_avg_blocklength,numel(obj.d_constellation)); m_reg = 0; - if obj.eq_avg_blocklength > 0 - averaging_window = zeros(obj.eq_avg_blocklength,1); - end - for k = 1:obj.sps:obj.x_length dc_cnt = dc_cnt+1; m=m+1; %get Sigal input vectors with correct length for VNLE x = obj.x_in(obj.Ne(1)+k-1:-1:k).'; -% - if obj.eq_avg_blocklength > 0 %% Das läuft gut mit 400er Fenster!! - averaging_window = circshift(averaging_window,obj.sps); - averaging_window(1:obj.sps,1) = x(1:obj.sps); - avg_(k) = mean(averaging_window); - x = x-avg_(k); - end - + %bring this signal to "special" VNLE format x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); @@ -270,61 +256,42 @@ classdef EQ_silas < handle x_d = [x_vnle;-d_vnle]; %Apply filter - if obj.mu_dc_dd > 0 - y(m) = obj.e_dc(end) + x_d.'* coeff; - else - y(m) = x_d.'* coeff; - end - - + y(m) = x_d.'* coeff; + %Decision 1 [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point - d_hat(k) = obj.d_constellation(symbol_idx); - - %Error between FFE & DFE filtered signal and Decision - obj.error(k) = y(m) - d_hat(k); -% - lvl_err(k,symbol_idx) = obj.error(k); + d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx); - lvl_err_mov(:,symbol_idx) = circshift(lvl_err_mov(:,symbol_idx),1); - lvl_err_mov(1,symbol_idx) = obj.error(k); + y_1(m,symbol_idx) = y(m); % after 1st iteration - %Decision 2 - y(m) = y(m)-mean(lvl_err_mov(:,symbol_idx),'omitnan'); - [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point - d_hat(k) = obj.d_constellation(symbol_idx); + %1st Error between FFE & DFE filtered signal and Decision + obj.error(m) = y(m) - d_hat(m,symbol_idx); + % lvl_err_1(m,symbol_idx) = y(m) - obj.d(m+1); + % + % %write current error to buffer + % lvl_err_mov(:,symbol_idx) = circshift(lvl_err_mov(:,symbol_idx),1); + % lvl_err_mov(1,symbol_idx) = obj.error(m); + % + % %Subtract a weighted error from y -> then Decision 2 + % err = mean(lvl_err_mov(:,symbol_idx),'omitnan'); + % + % y(m) = y(m)-(obj.mu_dc_dd(symbol_idx)*err); + % + % subtracted_error(m,symbol_idx) = obj.mu_dc_dd(symbol_idx)*mean(lvl_err_mov(:,symbol_idx),'omitnan'); + % + % y_2(m,symbol_idx) = y(m); + % + % [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision 2 for closest constellation point + % + % d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx); + % + % obj.error(m) = y(m) - d_hat(m,symbol_idx); + % + % lvl_err_2(m,symbol_idx) = y(m) - obj.d(m+1); %Update FFE and DFE coefficients - coeff = coeff - (mu_mat * (obj.error(k) * conj(x_d))); - - %Update DC error - dc_block(dc_cnt) = obj.error(k) ; - - if dc_cnt == obj.eq_parallelization_blocklength - if obj.eq_updatelatency > 1 - - obj.e_dc = circshift(obj.e_dc,1); - -% m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block))); -% -% obj.e_dc(1) = obj.e_dc(2) - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd); - obj.e_dc(1) = obj.e_dc(2) - sum(dc_block) .* obj.mu_dc_dd; - - else - %m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block))); - -% obj.e_dc = obj.e_dc - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd); - - obj.e_dc = obj.e_dc - sum(dc_block) .* obj.mu_dc_dd; - -% obj.e_dc = obj.e_dc - obj.mu_dc_dd * obj.error(k); %newapril - end - - dc_cnt = 0; - - end - + coeff = coeff - (mu_mat * (obj.error(m) * conj(x_d))); % Append new decision to decision feedback if obj.Nb(1) > 0 @@ -332,43 +299,52 @@ classdef EQ_silas < handle %shift up one index d_feedback(2:end) = d_feedback(1:end-1); %replace 1st index with current estimation - d_feedback(1) = d_hat(k); + d_feedback(1) = d_hat(m,symbol_idx); %build memorylike VNLE version d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); end - end end - %% -% -% b = movmean(lvl_err,[500 500],1,"omitnan"); -% -% figure(11) -% for i = 1:4 -% hold on -% stem(lvl_err(:,i)) -% end - %% - obj.y_out = (circshift( y.' ,-(obj.delay))).'; obj.d_out = d_hat(1:2:end); -% err = obj.error(1:2:end); -% res = NaN(8,length(err)); -% for lvl = 1:8 -% a = find(obj.d_out==obj.d_constellation(lvl)); -% res(lvl,a) = err(a); -% end -% mean(res,2,"omitnan"); -% -% figure(12) -% scatter(1:length(obj.y_out),obj.y_out,1,'.') -% hold on -% scatter(1:length(obj.d_out),obj.d_out,1,'.') + + % evm1 = mean(lvl_err_1,'omitnan'); + % + % evm2 = mean(lvl_err_2,'omitnan'); + % + % figure(112) + % stem(evm1,'LineStyle','--','Marker','square','LineWidth',1); + % hold on; + % stem(evm2,'LineStyle',':','Marker','v','LineWidth',1); + + + % figure(14) + % scatter(1:length(lvl_err_1),subtracted_error,1,'.') + % + % lvl_err___ = lvl_err_true(~isnan(lvl_err_true)); + % %lvl_err___ = lvl_err___-mean(lvl_err___); + % coeffs = arburg(lvl_err___,1000); + % fs_in = 92e9; + % [h,w] = freqz(1,coeffs,length(lvl_err___),"whole",fs_in); + % h = fftshift(h./max(abs(h))); + % freq_vec = linspace(-fs_in/2,fs_in/2,length(h)); + % figure(111) + % hold on + % plot(freq_vec.*1e-9,20*log10(h),'DisplayName','burg'); + + % + % spectrum_plot(y,92e9); + % + % d = 2^nextpow2(length(y)/16); + % + % figure(1111) + % hold on + % pwelch(y,hamming(d),d/2,d,92e9,"centered","power"); end diff --git a/Classes/04_DSP/EQ_silas_ofc.m b/Classes/04_DSP/EQ_silas_ofc.m new file mode 100644 index 0000000..48aa77d --- /dev/null +++ b/Classes/04_DSP/EQ_silas_ofc.m @@ -0,0 +1,456 @@ +classdef EQ_silas_ofc < handle + %EQ_SILAS FFE and DFE Equalizer Playground + + properties + % Important Signals + x_in %Input Sequence to be equalized + x_length + x_norm + + d %reference signal + d_norm + d_constellation %constellation points of the reference + + y_out %equalizer output signal + + % FFE coefficients always named with "e" + Ne + Ce %memory length FFE + Ie1 %Indice Combination of 1nd order FFE + Ie2 %Indice Combination of 2nd order FFE + Ie3 %Indice Combination of 3nd order FFE + e %coefficients for FFE + + % DFE coefficients always named with "b" + Nb + Cb %memory length DFE + Ib1 %Indice Combination of 1nd order DFE + Ib2 %Indice Combination of 2nd order DFE + Ib3 %Indice Combination of 3nd order DFE + b %coefficients for DFE + + error + e_ffe + e_dfe + e_dc + + % coefficients + mu_dc_train + mu_ffe_train + mu_dfe_train + + mu_dc_dd + mu_ffe_dd + mu_dfe_dd + mu_combined_dd % [1st order FFE, 2nd order FFE, 3rd order FFE, all orders DFE] + + delay + trainlength + sps + + trainloops + ddloops + + eq_parallelization_blocklength % block lengt of EQ (until now, only the dc subtraction is affected by this) + eq_updatelatency % time in symbols until the calculated updates reach the signal again (until now, only the dc subtraction is affected by this) + eq_avg_blocklength + + + end + + methods + function obj = EQ_silas_ofc(options) + %EQ_SILAS Construct an instance of this class + arguments(Input) + options.Ne = [50 5 0] %Number of FFE coefficients (1st, 2nd and 3rd order) + options.Nb = [30 5 3] %Number of DFE coefficients (1st, 2nd and 3rd order) + options.trainloops = 2; + options.trainlength = 4096; + options.ddloops = 2; + + options.delay = 0; + options.sps = 2; + + options.mu_dc_train = 0.01; + options.mu_ffe_train = 0.005; + options.mu_dfe_train = 0.005; + + options.mu_dc_dd = 0.01; + options.mu_ffe_dd = [0.0004 0.0005 0.0006]; + options.mu_dfe_dd = 0.0005; + + options.eq_parallelization_blocklength = 1; + options.eq_updatelatency = 1; + options.eq_avg_blocklength = 0; + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + % Generate helpful vectors and initialize the filters with + % correct length: + + obj.Ce = obj.calcVNLEMemoryLength(obj.Ne); + + [obj.Ie2,obj.Ie3] = obj.calcIndiceVectors(obj.Ne); + + obj.e = zeros(sum(obj.Ce),1); + + + obj.Cb = obj.calcVNLEMemoryLength(obj.Nb); + + [obj.Ib2,obj.Ib3] = obj.calcIndiceVectors(obj.Nb); + + obj.b = zeros(sum(obj.Cb),1); + + end + + function [signalclass_out] = process(obj,signalclass_in, reference_signalclass_in) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + signalclass_in = signalclass_in.normalize("mode","rms"); + + % Process the EQ optimization + obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); + + signalclass_in.signal = obj.y_out'; + % append to logbook + lbdesc = ['EQ von Silas ist gelaufen ']; + signalclass_in = signalclass_in.logbookentry(lbdesc); + + % write to output + signalclass_out = signalclass_in; + + end + + function process_(obj,x_in,d_in) + + % 1) prepare signals + obj.e_dc = mean(x_in); + + % 1.1) Input Signal + obj.x_in = [zeros(1,floor(obj.Ne(1)/2)) x_in zeros(1,obj.Ne(1))]; + obj.x_length = length(x_in); + obj.x_norm = obj.calcPowerNormalization(x_in); + + % 1.2 Reference Signal // Constellation + obj.d = [zeros(1,obj.Nb(1)-1) d_in zeros(1,obj.Nb(1))]; + obj.d_constellation = unique(d_in); + obj.d_norm = obj.calcPowerNormalization(d_in); + + % 1.3 Training + obj.trainingMode(); + + % 1.4 Decision Directed Mode + obj.decisionDirectedMode(); + + end + + %% Adaptive Equalization Modes + + function trainingMode(obj) + + dc_block = ones(obj.eq_parallelization_blocklength,1); + + for tloop = 1:obj.trainloops + m = 1+obj.delay; + dc_cnt = 0; + for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength + m = m+1; + dc_cnt = dc_cnt+1; + + %get Sigal input vectors with correct length for VNLE + x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).'; + + x_in_vnle_format = obj.calcVNLENonlinVecs(x_in_block,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); + + %get Reference input vectors with correct length for VNLE + d_block = obj.d(obj.Nb(1)-obj.delay+m-2:-1:m-obj.delay-1).'; + d_vnle_format = obj.calcVNLENonlinVecs(d_block,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + + obj.e_ffe = obj.e.' * x_in_vnle_format; + + obj.e_dfe = obj.b.' * d_vnle_format; + + % Calculate the Error + obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); + + if obj.mu_ffe_train ~= 0 + %update FFE coefficients with LMS + obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train; + else + %update FFE coefficients with NLMS + obj.e = obj.e - obj.error*x_in_vnle_format/(x_in_vnle_format.'*x_in_vnle_format); + end + + %update DFE coefficients with LMS + obj.b = obj.b + obj.mu_dfe_train*obj.error*d_vnle_format; + + %update DC error + dc_block(dc_cnt) = obj.error .* obj.mu_dc_train; + + if dc_cnt == obj.eq_parallelization_blocklength + obj.e_dc = obj.e_dc - mean(dc_block(dc_cnt)); + dc_cnt = 0; + end + + end + + end + + end + + + function decisionDirectedMode(obj) + + %start the dd mode with coefficients from training + coeff = [obj.e;obj.b]; + obj.e_dc = ones(obj.eq_updatelatency,1).*obj.e_dc; + dc_block = ones(obj.eq_parallelization_blocklength,1); + lvl_err_true = NaN(length(obj.d),numel(obj.d_constellation)); + for ddloop = 1:obj.ddloops + + m = 0; + dc_cnt = 0; + + mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe + ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe + ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)... %3rd order ffe + ones(1,sum(obj.Cb))*obj.mu_dfe_dd]); %all order dfe + + mu_ffe = [ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe + ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe + ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)]; + + mu_dfe = ones(1,sum(obj.Cb))*obj.mu_dfe_dd; + + y = zeros(1,floor(obj.x_length/obj.sps)); + d_feedback = zeros(obj.Cb(1),1); + d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + d_hat = zeros(obj.x_length,1); + + m_reg = 0; + + if obj.eq_avg_blocklength > 0 + averaging_window = zeros(obj.eq_avg_blocklength,1); + end + + for k = 1:obj.sps:obj.x_length + dc_cnt = dc_cnt+1; + m=m+1; + + %get Sigal input vectors with correct length for VNLE + x = obj.x_in(obj.Ne(1)+k-1:-1:k).'; +% + if obj.eq_avg_blocklength > 0 %% Das läuft gut mit 400er Fenster!! + averaging_window = circshift(averaging_window,obj.sps); + averaging_window(1:obj.sps,1) = x(1:obj.sps); + avg_(k) = mean(averaging_window); + x = x-avg_(k); + end + + %bring this signal to "special" VNLE format + x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm); + + %combine FFE with DFE to one vector (cursor between the two sequences) + x_d = [x_vnle;-d_vnle]; + + + %Apply filter + %y(m) = (m_reg(end)*dc_cnt + obj.e_dc(end)) + x_d.'* coeff; + if obj.mu_dc_dd > 0 + y(m) = obj.e_dc(end) + x_d.'* coeff; + else + y(m) = x_d.'* coeff; + end + +% if obj.eq_avg_blocklength > 0 %% Das läuft nicht gut!! +% averaging_window = circshift(averaging_window,obj.sps); +% averaging_window(1:obj.sps,1) = y(m); +% avg_(m) = mean(averaging_window); +% y(m) = y(m)-avg_(m); +% end + + %Decision + [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point + d_hat(k) = obj.d_constellation(symbol_idx); + + %Error between FFE & DFE filtered signal and Decision + obj.error(k) = y(m) - d_hat(k); + lvl_err_true(m,symbol_idx) = y(m) - obj.d(m+1); + +% if obj.eq_avg_blocklength > 0 %% Das läuft nicht gut!! +% averaging_window = circshift(averaging_window,obj.sps); +% averaging_window(1:obj.sps,1) = y(m); +% avg_(m) = mean(averaging_window); +% y(m) = y(m)-avg_(m); +% end + + %Update FFE and DFE coefficients + coeff = coeff - mu_mat*obj.error(k) * conj(x_d); + + %Update DC error + dc_block(dc_cnt) = obj.error(k) ; + + if dc_cnt == obj.eq_parallelization_blocklength + if obj.eq_updatelatency > 1 + + obj.e_dc = circshift(obj.e_dc,1); + +% m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block))); +% +% obj.e_dc(1) = obj.e_dc(2) - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd); + obj.e_dc(1) = obj.e_dc(2) - sum(dc_block) .* obj.mu_dc_dd; + + else + %m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block))); + +% obj.e_dc = obj.e_dc - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd); + + obj.e_dc = obj.e_dc - sum(dc_block) .* obj.mu_dc_dd; + end + + dc_cnt = 0; + + end + + + % Append new decision to decision feedback + if obj.Nb(1) > 0 + + %shift up one index + d_feedback(2:end) = d_feedback(1:end-1); + %replace 1st index with current estimation + d_feedback(1) = d_hat(k); + %build memorylike VNLE version + d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm); + + end + + end + end + + obj.y_out = (circshift( y.' ,-(obj.delay))).'; + + end + + %% Functions needed During Adaption + function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) + % These are the second and third order input signal products of the VNLE EQ + % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) + + x1 = x_in_block; + x2 = []; + x3 = []; + + if N_(2) > 0 + delta_2 = round((N_(1)-N_(2))/2); + input_vec_se = x_in_block(delta_2:end)/norm_(2); %TODO normalization step + x2 = input_vec_se(I_2(:,1)).*input_vec_se(I_2(:,2)); + end + + if N_(3) > 0 + delta_3 = round((N_(1)-N_(3))/2); + input_vec_th = x_in_block(delta_3:end)/norm_(3); + x3 = input_vec_th(I_3(:,1)).*input_vec_th(I_3(:,2)).*input_vec_th(I_3(:,3)); + end + + x_in_vnle_format = [x1;x2;x3]; + + end + + + %% Functions needed for Preparation + function [C] = calcVNLEMemoryLength(~,N) + + %calculates the memory length of VNLE + C = zeros(size(N)); + + for o = 1:numel(N) + switch o + case 1 + C(o) = N(o); + case 2 + C(o) = N(o)*(N(o)+1) / 2; + case 3 + C(o) = N(o)*(N(o)+1)*(N(o)+2) / 6; + end + end + + end + + function [indvec2nd, indvec3rd] = calcIndiceVectors(~,N) + + % Init vectors of 2nd and 3rd order coefficient indices -> + % yield combination with + + for order = 2:numel(N) + n = N(order); + v = 1:n; % Ursprünglicher Vektor + row = 1; + + % Schleifen zur Generierung des Indize Vektors + switch order + + case 2 + + indvec2nd = zeros(n*(n+1)/2, order); + for i = 1:n + for j = i:n + indvec2nd(row, :) = [v(i) v(j)]; + row = row + 1; + end + end + + case 3 + + indvec3rd = zeros(n*(n+1)*(n+2)/6, 3); + for i = 1:n + for j = i:n + for k = j:n + indvec3rd(row, :) = [v(i) v(j) v(k)]; + row = row + 1; + end + end + end + end + end + + end + + function powerNorm = calcPowerNormalization(~,v) + + powerNorm(1) = sqrt(mean(abs(v ).^2)); + powerNorm(2) = sqrt(mean(abs(v.^2).^2)); + powerNorm(3) = sqrt(mean(abs(v.^3).^2)); + + end + + + end +end + + + + + + + + + + + + + + + + + + + + + + diff --git a/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m b/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m index cf1bde1..3496315 100644 --- a/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m +++ b/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m @@ -54,8 +54,6 @@ classdef EQ_silas_sliding_window_dc_removal < handle eq_blocklength % block lengt of EQ (until now, only the dc subtraction is affected by this) eq_updatelatency % time in symbols until the calculated updates reach the signal again (until now, only the dc subtraction is affected by this) - - end methods @@ -110,7 +108,7 @@ classdef EQ_silas_sliding_window_dc_removal < handle % actual processing of the signal (steps 1. - 3.) % 1 normalize RMS - %signalclass_in = signalclass_in.normalize("mode","rms"); + signalclass_in = signalclass_in.normalize("mode","rms"); % Process the EQ optimization obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); diff --git a/Classes/04_DSP/FFE.m b/Classes/04_DSP/FFE.m new file mode 100644 index 0000000..9b39be8 --- /dev/null +++ b/Classes/04_DSP/FFE.m @@ -0,0 +1,170 @@ +classdef FFE < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + properties + sps % usually 2 + order + e + error + + len_tr + mu_tr + epochs_tr + + mu_dd + epochs_dd + + constellation + + decide + end + + methods + function obj = FFE(options) + arguments(Input) + + options.sps = 2; + options.order = 15; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + options.decide = false; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.order,1); + obj.error = 0; + + end + + function [X] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + + % Decision Directed Mode + N = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + + end + + function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz) + + arguments + obj + x + d + mio + epochs + N + training + showviz + end + + x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + + if showviz + f = figure(111); + subplot(2,2,1:2); + hold on + a = scatter(1:numel(x),x,1,'.'); + a2 = scatter(1,1,1,'.'); + a3 = scatter(1,1,2,'.'); + a4 = xline(1); + ylim([-3 3]) + xlim([0 length(x)]); + % subplot(2,2,3) + % dplot = x(1:1+500); + % b = scatter(1:numel(dplot),dplot,5,'x'); + % xline(1) + % xline(obj.order) + % ylim([-3 3]) + % xlim([0 500]); + subplot(2,2,3:4) + c = stem(obj.e); + ylim([-1 1]) + drawnow + end + + for epoch = 1 : epochs + symbol = 0; + for sample = 1 : obj.sps : N + + symbol = symbol+1; + + U = x(obj.order+sample-1:-1:sample); + + y(symbol,1) = obj.e.' * U; % Calculating output of LMS __ * | + + if training + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + + if mio ~= 0 + obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS + else + normalizationfactor = (U.' * U); + obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS + end + + if mod(sample,100) == 1 && showviz + a2.XData = 1:2*numel(y); + a2.YData = repelem(y, 2); + a3.XData = 1:2*numel(d_hat); + a3.YData = repelem(d_hat, 2); + a4.Value = sample; + % b.YData = x(symbol:symbol+500); + c.YData = obj.e; + drawnow; + end + + obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + + end + end + + end + + end +end + diff --git a/Classes/04_DSP/FFE_DFE.m b/Classes/04_DSP/FFE_DFE.m new file mode 100644 index 0000000..0841f8d --- /dev/null +++ b/Classes/04_DSP/FFE_DFE.m @@ -0,0 +1,194 @@ +classdef FFE_DFE < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + % Eq = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",21,"dfe_order",0,"sps",2,"decide",1); + + properties + sps % usually 2 + ffe_order + dfe_order + e + b + error + + len_tr + ffe_mu_tr + dfe_mu_tr + epochs_tr + + ffe_mu_dd + dfe_mu_dd + epochs_dd + + constellation + + decide + end + + methods + function obj = FFE_DFE(options) + arguments(Input) + + options.sps = 2; + + options.ffe_order = 15; + options.dfe_order = 2; + + options.len_tr = 4096; + options.ffe_mu_tr = 0; + options.dfe_mu_tr = 0; + options.epochs_tr = 5; + + options.ffe_mu_dd = 1e-5; + options.dfe_mu_dd = 1e-5; + options.epochs_dd = 5; + + options.decide = false; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + obj.e = zeros(obj.ffe_order,1); + obj.b = zeros(obj.dfe_order,1); + obj.error = 0; + + end + + function [X] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.ffe_mu_tr,obj.dfe_mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + + % Decision Directed Mode + N = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.ffe_mu_dd,obj.dfe_mu_dd,obj.epochs_dd,N,training,showviz); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.ffe_order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + + end + + function [y,d_hat] = equalize(obj,x,d,ffe_mu,dfe_mu,epochs,N,training,showviz) + + arguments + obj + x + d + ffe_mu + dfe_mu + epochs + N + training + showviz + end + + + mu = diag([ones(1,obj.ffe_order(1))*ffe_mu(1) ... + ones(1,obj.dfe_order(1))*dfe_mu(1) ]); + + + x = [zeros(floor(obj.ffe_order/2),1); x; zeros(obj.ffe_order,1)]; + %d = [zeros(obj.dfe_order-1,1); d; zeros(obj.dfe_order,1)]; + d_ = zeros(obj.dfe_order(1),1); + coeff = [obj.e;obj.b]; + + if showviz + f = figure(111); + subplot(2,2,1:2); + hold on + a = scatter(1:numel(x),x,1,'.'); + a2 = scatter(1,1,1,'.'); + a3 = scatter(1,1,2,'.'); + a4 = xline(1); + ylim([-3 3]) + xlim([0 length(x)]); + subplot(2,2,3:4) + c = stem(obj.e); + ylim([-1 1]) + drawnow + end + + for epoch = 1 : epochs + symbol = 0; + for sample = 1 : obj.sps : N + + symbol = symbol+1; + + x_ = x(obj.ffe_order+sample-1:-1:sample); + v = [x_;d_]; + + y(symbol,1) = coeff.' * v; % Calculating output of LMS __ * | + + if training + d_hat(symbol,1) = d(symbol); + else + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + end + + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + + if ~all(mu == 0,'all') %not all mu values are zero + coeff = coeff - (mu * err(symbol) * v) ; % Weight update rule of LMS + else + normalizationfactor = (v.' * v); + coeff = coeff - err(symbol) * v / normalizationfactor; % Weight update rule of NLMS + end + + % Append new decision to decision feedback + if obj.dfe_order(1) > 0 + %shift up one index + d_(2:end) = d_(1:end-1); + %replace 1st index with current estimation + d_(1) = d_hat(symbol); + end + + if mod(sample,100) == 1 && showviz + a2.XData = 1:2*numel(y); + a2.YData = repelem(y, 2); + a3.XData = 1:2*numel(d_hat); + a3.YData = repelem(d_hat, 2); + a4.Value = sample; + c.YData = obj.e; + drawnow; + end + + obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + + end + end + + obj.e = coeff(1:obj.ffe_order); + obj.b = coeff(obj.ffe_order+1:end); + + + + end + + end +end + diff --git a/Classes/04_DSP/VNLE.m b/Classes/04_DSP/VNLE.m new file mode 100644 index 0000000..60dc09d --- /dev/null +++ b/Classes/04_DSP/VNLE.m @@ -0,0 +1,293 @@ +classdef VNLE < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + properties + sps % usually 2 + order + e + error + + len_tr + mu_tr + epochs_tr + + mu_dd + epochs_dd + + constellation + + decide + + x_norm + ce + ie2 + ie3 + end + + methods + function obj = VNLE(options) + arguments(Input) + + options.sps = 2; + options.order = [15,2,2]; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.mu_dd = 1e-5; + options.epochs_dd = 5; + + options.decide = false; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + + obj.error = 0; + + end + + function [X] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + obj.x_norm = obj.calcPowerNormalization(X.signal); + obj.ce = obj.calcVNLEMemoryLength(obj.order); + [obj.ie2,obj.ie3] = obj.calcIndiceVectors(obj.order); + + obj.e = zeros( sum(obj.ce) ,1); + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + + % Decision Directed Mode + N = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + + end + + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) + + arguments + obj + x + d + mu + epochs + N + training + showviz + end + + if all(mu == mu(1)) + % mu = mu(1); + mu = diag(ones(1,sum(obj.ce))*mu(1)); + else + mu = diag([ones(1,obj.ce(1))*mu(1) ... + ones(1,obj.ce(2))*mu(2) ... + ones(1,obj.ce(3))*mu(3) ]); + end + + x = [zeros(floor(obj.order(1)/2),1); x; zeros(obj.order(1),1)]; + + if showviz + f = figure(111); + subplot(2,2,1:2); + hold on + a = scatter(1:numel(x),x,1,'.'); + a2 = scatter(1,1,1,'.'); + a3 = scatter(1,1,2,'.'); + a4 = xline(1); + ylim([-3 3]) + xlim([0 length(x)]); + subplot(2,2,3:4) + c = stem(obj.e); + ylim([-1 1]) + drawnow + end + + for epoch = 1 : epochs + symbol = 0; + for sample = 1 : obj.sps : N + + symbol = symbol+1; + + + % x_in = x(obj.order(1)+sample+(obj.sps-1):-1:sample+obj.sps); + x_in = x(obj.order(1)+sample-1:-1:sample); + x_in = obj.calcVNLENonlinVecs(x_in,obj.ie2,obj.ie3,obj.order,obj.x_norm); + + y(symbol,1) = obj.e.' * x_in; % Calculating output of LMS __ * | + + if training + err(symbol) = y(symbol) - d(symbol); % Instantaneous error + else + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + end + + if ~all(mu==0,'all') %mu has not only zeros + obj.e = obj.e - (mu * err(symbol) * x_in) ; % Weight update rule of LMS + else + normalizationfactor = (x_in.' * x_in); + obj.e = obj.e - err(symbol) * x_in / normalizationfactor; % Weight update rule of NLMS + end + + if mod(sample,100) == 1 && showviz + a2.XData = 1:2*numel(y); + a2.YData = repelem(y, 2); + a3.XData = 1:2*numel(d_hat); + a3.YData = repelem(d_hat, 2); + a4.Value = sample; + % b.YData = x(symbol:symbol+500); + c.YData = obj.e; + drawnow; + end + + obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + + end + end + + end + + %% Functions needed During Adaption + function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) + % These are the second and third order input signal products of the VNLE EQ + % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) + l1=length(x_in_block); + l2=length(I_2); + l3=length(I_3); + final_length = l1+l2+l3; + + x_in_vnle_format = zeros(final_length,1); + + idx = l1; + x_in_vnle_format(1:idx) = x_in_block; + + if N_(2) > 0 + delta_2 = round((N_(1)-N_(2)) / 2); + input_vec_se = x_in_block(delta_2:end) / norm_(2); %TODO normalization step + + % Extract columns from I_2 + col1 = input_vec_se(I_2(:,1)); + col2 = input_vec_se(I_2(:,2)); + + x2 = col1 .* col2; + x_in_vnle_format(idx+1:idx+l2) = x2; + end + + if N_(3) > 0 + delta_3 = round((N_(1)-N_(3))/2); + input_vec_th = x_in_block(delta_3:end) / norm_(3); + + % Extract columns from I_3 + col1 = input_vec_th(I_3(:,1)); + col2 = input_vec_th(I_3(:,2)); + col3 = input_vec_th(I_3(:,3)); + + % Perform matrix multiplication + x3 = col1 .* col2 .* col3; + + idx = idx+l2; + x_in_vnle_format(idx+1:idx+l3) = x3; + end + + end + + %% Functions needed for Preparation + function [C] = calcVNLEMemoryLength(~,N) + + %calculates the memory length of VNLE + C = zeros(size(N)); + + for o = 1:numel(N) + switch o + case 1 + C(o) = N(o); + case 2 + C(o) = N(o)*(N(o)+1) / 2; + case 3 + C(o) = N(o)*(N(o)+1)*(N(o)+2) / 6; + end + end + + end + + function [indvec2nd, indvec3rd] = calcIndiceVectors(~,N) + + % Init vectors of 2nd and 3rd order coefficient indices -> + % yield combination with + indvec2nd=[]; + indvec3rd=[]; + for o = 2:numel(N) + n = N(o); + v = 1:n; % Ursprünglicher Vektor + row = 1; + + % Schleifen zur Generierung des Indize Vektors + switch o + + case 2 + + indvec2nd = zeros(n*(n+1)/2, o); + for i = 1:n + for j = i:n + indvec2nd(row, :) = [v(i) v(j)]; + row = row + 1; + end + end + + case 3 + + indvec3rd = zeros(n*(n+1)*(n+2)/6, 3); + for i = 1:n + for j = i:n + for k = j:n + indvec3rd(row, :) = [v(i) v(j) v(k)]; + row = row + 1; + end + end + end + end + end + + end + + function powerNorm = calcPowerNormalization(~,v) + + powerNorm(1) = sqrt(mean(abs(v ).^2)); + powerNorm(2) = sqrt(mean(abs(v.^2).^2)); + powerNorm(3) = sqrt(mean(abs(v.^3).^2)); + + end + + end +end + diff --git a/Classes/matlab.mat b/Classes/matlab.mat new file mode 100644 index 0000000..e33da65 Binary files /dev/null and b/Classes/matlab.mat differ diff --git a/Datatypes/signalform.m b/Datatypes/signalform.m new file mode 100644 index 0000000..b6a325b --- /dev/null +++ b/Datatypes/signalform.m @@ -0,0 +1,10 @@ +classdef signalform < int32 + + enumeration + sine (1) + sawtooth (2) + square (3) + noise (4) + end + +end \ No newline at end of file diff --git a/projects/AWG_characterization/output_power.m b/projects/AWG_characterization/output_power.m new file mode 100644 index 0000000..2fb812c --- /dev/null +++ b/projects/AWG_characterization/output_power.m @@ -0,0 +1,195 @@ + +M=4; + +fdac = 256e9;%fsym; +fadc = 256e9; +fsym = ([32:16:240].*1e9); + +%fsym = 160e9; +% 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)); +end + +% 2 ) Build Inf. signal class +bits = Informationsignal(bitpattern); + +% 5) AWG (lowpass, quantization, sample and hold) +kover = 8; +LP_awg = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*kover,"filterType",filtertypes.butterworth); + +powerlist = []; +input_papr = zeros(numel(fsym),1); +output_papr = zeros(numel(fsym),1); + +for i = 1:length(fsym) + + %%%%% Map to PAM %%%%%% + digimod_out = PAMmapper(M,0).map(bits); + digimod_out.fs = fsym(i); + + % Y.plot("fignum",2,'displayname',['PAM 4 Signal; Raised Cosine Alpha: ',num2str(rrca)]); + % X.plot("fignum",2,'displayname',['PAM 4 Signal; Raised Cosine Alpha: ',num2str(rrca)]) + % digimod_out.plot("fignum",2,'displayname',['PAM 4 Signal']) + + %X.spectrum("displayname",['PAM 4; Baudrate: ',num2str(fsym(i).*1e-9), ' GBd', num2str(rrca)],'fignum',4); + %%%%% AWG %%%%%% + + X = digimod_out; + + X1 = M8199B("kover",kover).process(X); + + X2 = M8199A("kover",kover).process(X); + + X3 = M8196A("kover",kover).process(X); + + % X.spectrum("displayname",['Baudrate: ',num2str(fsym(i).*1e-9), ' GBd'],'fignum',4); + + powerlist1(i) = X1.power; + vpplist1(i) = max(X1.signal)-min(X1.signal); + paprlist1(i) = X1.papr_lin; + + powerlist2(i) = X2.power; + vpplist2(i) = max(X2.signal)-min(X2.signal); + paprlist2(i) = X2.papr_lin; + + if fsym(i) <= 113e9 + powerlist3(i) = X3.power; + vpplist3(i) = max(X3.signal)-min(X3.signal); + paprlist3(i) = X3.papr_lin; + else + powerlist3(i) =NaN; + vpplist3(i) = NaN; + paprlist3(i) = NaN; + end + + + + %%%%% Pulseforming %%%%%% + rrca=0.3; + X = Pulseformer("fsym",fsym(i),"fdac",256e9,"pulse","rrc","pulselength",16,"rrcalpha",rrca).process(digimod_out); + + % % %%%%% Clip to PAM range %%%%%% + min_ = min(digimod_out.signal).*1.3; + max_ = max(digimod_out.signal).*1.3; + X.signal = clip(X.signal,min_,max_); + + X1 = M8199B("kover",kover).process(X); + + X2 = M8199A("kover",kover).process(X); + + X = Pulseformer("fsym",fsym(i),"fdac",92e9,"pulse","rrc","pulselength",16,"rrcalpha",rrca).process(digimod_out); + + % % %%%%% Clip to PAM range %%%%%% + min_ = min(digimod_out.signal).*1.3; + max_ = max(digimod_out.signal).*1.3; + X.signal = clip(X.signal,min_,max_); + + X3 = M8196A("kover",kover).process(X); + + powerlist1_opt(i) = X1.power; + vpplist1_opt(i) = max(X1.signal)-min(X1.signal); + paprlist1_opt(i) = X1.papr_lin; + + powerlist2_opt(i) = X2.power; + vpplist2_opt(i) = max(X2.signal)-min(X2.signal); + paprlist2_opt(i) = X2.papr_lin; + + if fsym(i) <= 113e9 + powerlist3_opt(i) = X3.power; + vpplist3_opt(i) = max(X3.signal)-min(X3.signal); + paprlist3_opt(i) = X3.papr_lin; + else + powerlist3_opt(i) =NaN; + vpplist3_opt(i) = NaN; + paprlist3_opt(i) = NaN; + end + +end + +cols = linspecer(3); + +figure(6) +hold on +scatter(fsym.*1e-9,vpplist1,'DisplayName','M8199B','MarkerFaceColor',cols(1,:),'MarkerEdgeColor',cols(1,:),'LineWidth',2); +scatter(fsym.*1e-9,vpplist2,'DisplayName','M8199A','MarkerFaceColor',cols(2,:),'MarkerEdgeColor',cols(2,:),'LineWidth',2); +scatter(fsym.*1e-9,vpplist3,'DisplayName','M8196A','MarkerFaceColor',cols(3,:),'MarkerEdgeColor',cols(3,:),'LineWidth',2); +xticks(fsym.*1e-9) +xlabel("Baudrate in GBaud"); +ylabel("Vpp in V") +legend +ylim([0.3 1.4]) +thickenfigure; + +figure(7) +hold on +scatter(fsym.*1e-9,powerlist1,'DisplayName','M8199B','MarkerFaceColor',cols(1,:),'MarkerEdgeColor',cols(1,:),'LineWidth',2); +scatter(fsym.*1e-9,powerlist2,'DisplayName','M8199A','MarkerFaceColor',cols(2,:),'MarkerEdgeColor',cols(2,:),'LineWidth',2); +scatter(fsym.*1e-9,powerlist3,'DisplayName','M8196A','MarkerFaceColor',cols(3,:),'MarkerEdgeColor',cols(3,:),'LineWidth',2); +xticks(fsym.*1e-9) +xlabel("Baudrate in GBaud"); +ylabel("Output Power in dBm") +legend +ylim([-12 4]) +thickenfigure; + +figure(8) +hold on +scatter(fsym.*1e-9,paprlist1,'DisplayName','M8199B','MarkerFaceColor',cols(1,:),'MarkerEdgeColor',cols(1,:),'LineWidth',2); +scatter(fsym.*1e-9,paprlist2,'DisplayName','M8199A','MarkerFaceColor',cols(2,:),'MarkerEdgeColor',cols(2,:),'LineWidth',2); +scatter(fsym.*1e-9,paprlist3,'DisplayName','M8196A','MarkerFaceColor',cols(3,:),'MarkerEdgeColor',cols(3,:),'LineWidth',2); +xticks(fsym.*1e-9) +xlabel("Baudrate in GBaud"); +ylabel("PAPR linear") +legend +ylim([3 10]) +thickenfigure; + + +figure(9) +hold on +scatter(fsym.*1e-9,vpplist1_opt,'DisplayName','M8199B','Marker','x','MarkerFaceColor',cols(1,:),'MarkerEdgeColor',cols(1,:),'LineWidth',3); +scatter(fsym.*1e-9,vpplist2_opt,'DisplayName','M8199A','Marker','x','MarkerFaceColor',cols(2,:),'MarkerEdgeColor',cols(2,:),'LineWidth',3); +scatter(fsym.*1e-9,vpplist3_opt,'DisplayName','M8196A','Marker','x','MarkerFaceColor',cols(3,:),'MarkerEdgeColor',cols(3,:),'LineWidth',3); +xticks(fsym.*1e-9) +xlabel("Baudrate in GBaud"); +ylabel("Vpp in V") +legend +ylim([0.3 1.4]) +thickenfigure; + +figure(10) +hold on +scatter(fsym.*1e-9,powerlist1_opt,'DisplayName','M8199B','Marker','x','MarkerFaceColor',cols(1,:),'MarkerEdgeColor',cols(1,:),'LineWidth',3); +scatter(fsym.*1e-9,powerlist2_opt,'DisplayName','M8199A','Marker','x','MarkerFaceColor',cols(2,:),'MarkerEdgeColor',cols(2,:),'LineWidth',3); +scatter(fsym.*1e-9,powerlist3_opt,'DisplayName','M8196A','Marker','x','MarkerFaceColor',cols(3,:),'MarkerEdgeColor',cols(3,:),'LineWidth',3); +xticks(fsym.*1e-9) +xlabel("Baudrate in GBaud"); +ylabel("Output Power in dBm") +legend +ylim([-12 4]) +thickenfigure; + +figure(11) +hold on +scatter(fsym.*1e-9,paprlist1_opt,'DisplayName','M8199B','Marker','x','MarkerFaceColor',cols(1,:),'MarkerEdgeColor',cols(1,:),'LineWidth',3); +scatter(fsym.*1e-9,paprlist2_opt,'DisplayName','M8199A','Marker','x','MarkerFaceColor',cols(2,:),'MarkerEdgeColor',cols(2,:),'LineWidth',3); +scatter(fsym.*1e-9,paprlist3_opt,'DisplayName','M8196A','Marker','x','MarkerFaceColor',cols(3,:),'MarkerEdgeColor',cols(3,:),'LineWidth',3); +xticks(fsym.*1e-9) +xlabel("Baudrate in GBaud"); +ylabel("PAPR linear") +legend +ylim([3 10]) +thickenfigure; + +autoArrangeFigures diff --git a/projects/MPI_April/auswertung_1.m b/projects/MPI_April/auswertung_1.m index 1c8275d..ff889e9 100644 --- a/projects/MPI_April/auswertung_1.m +++ b/projects/MPI_April/auswertung_1.m @@ -1,7 +1,4 @@ - - - vp = wh.parameter.vp.values(2); vb = wh.parameter.vb.values(1); rop = wh.parameter.rop.values; diff --git a/projects/MPI_April/mpi_simulation_cspr.m b/projects/MPI_April/mpi_simulation_cspr.m index 9dca01d..811b65b 100644 --- a/projects/MPI_April/mpi_simulation_cspr.m +++ b/projects/MPI_April/mpi_simulation_cspr.m @@ -1,17 +1,17 @@ %% Settings clear -for M=[4] - +for M=[8] filename = '112G_2'; load_sequence = 0; - datarate = 448e9; + datarate = 224e9; - kover = 8; + kover = 16; fsym = round(datarate*1e-9 / log2(M))*1e9; - fdac = fsym;%256e9; + + fdac = 256e9;%fsym; fadc = 256e9; lowpass_cutoff = fsym/2 * 1.1; @@ -20,11 +20,11 @@ for M=[4] phd_bw = lowpass_cutoff; scp_bw = lowpass_cutoff; - 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); + LP_awg = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true); + LP_modulator= Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true); + LP_opt = Filter('filtdegree',3,"f_cutoff",fsym/log2(M).*1.5,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + LP_phd = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true); + LP_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); % 1) PRBS Generation O = 18; %order of prbs @@ -48,15 +48,6 @@ for M=[4] digimod_out = PAMmapper(M,0).map(bits); digimod_out.fs = fsym; - - % linearGain = 1; - % limit = 1; - % SaturatingAmplifier = serdes.SaturatingAmplifier('Mode',1,... - % 'Limit',limit,'LinearGain',linearGain); - % X.signal = SaturatingAmplifier(X.signal); - % X.signal = min(max(X.signal,-0.8),0.8); - % X = X.normalize("mode","oneone"); - sir = [20:2:36]; %decibel = attenuation of interference path laser_linewidth = [1e5 1e6 10e6]; pn_key = [1:10]; @@ -64,12 +55,12 @@ for M=[4] vb = [1:0.1:1.8]; rop = 0; - sir = 28; - laser_linewidth = 10e6; + sir = 25; + laser_linewidth = 0; pn_key = 9; - vp = 0.5; + vp = 1;%0.5; vb = 1;%[1:0.1:1.8]; - mpi_path = 50; + mpi_path = 0; cnt = 1; @@ -85,62 +76,14 @@ for M=[4] %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 OPTICAL DOMAIN u_pi = 2; @@ -179,12 +122,12 @@ for M=[4] Combined_sig.signal = Combined_sig.signal(ceil(dly):end); % 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); + Combined_sig = Fiber("fsimu",Combined_sig.fs,"fiber_length",0,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.08).process(Combined_sig); 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); + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop(i)).process(Combined_sig); rop_save(s,l,pnk,n,m,i) = Rx_sig.power; % Square Law @@ -204,18 +147,8 @@ for M=[4] % Sync Rx signal with reference [Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",digimod_out,"fs_ref",fsym); - - % % % simple EQ (optimum mudc: 0.05 -> 0.005) - % EQ_sig = EQ_silas_plain("Ne",[20,8,8],"Nb",[2,0,0],"trainlength",4096,"mu_dc_dd",0.005,"mu_dc_train",0.05,... - % "mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0006 0.0003 0.005],"ddloops",3,'trainloops',3,'sps',2).process(Rx_sig,digimod_out); - -% EQ_sig = EQ("K",2,"plottrain",0,"plotfinal",0,... -% "training_length",4096,"training_loops",3,... -% "Ne",[50,8,8],"Nb",[2,0,0],... -% "DCmu",0.005,"DDmu",[0.0004 0.0006 0.0003 0.005],"DFEmu",0.005,"FFEmu",0.00,... -% "dd_loops",3,"epsilon",[10 100 1000 ],"M",2,... -% "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0).process(Scpe_sig,digimod_out); - + + Scpe_sig.spectrum; [EQ_sig,EQ_sym] = EQ_silas("Ne",[50,8,8],"Nb",[2,0,0],"trainlength",4096,... "sps",2,... "mu_dc_dd",0.00,... @@ -227,12 +160,12 @@ for M=[4] "ddloops",3,... "trainloops",3,... "eq_parallelization_blocklength",1, ... - "eq_updatelatency",1,... + "eq_updatelatency",0,... "eq_avg_blocklength",0).process(Scpe_sig,digimod_out); % Demap - Rx_Bits = PAMmapper(M,0).demap(EQ_sym); + Rx_Bits = PAMmapper(M,0).demap(EQ_sig); % 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); @@ -240,7 +173,7 @@ for M=[4] 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 diff --git a/projects/MPI_August/channel_model_mpi.m b/projects/MPI_August/channel_model_mpi.m new file mode 100644 index 0000000..19fd984 --- /dev/null +++ b/projects/MPI_August/channel_model_mpi.m @@ -0,0 +1,28 @@ +function opt_signal = channel_model_mpi(opt_signal,link_total_meter,oneway_interference_meter,sir) + +%%%%% Local Parameter %%%%%%%% + + + +%%%%% Ping Pong/ Interference Path %%%%%% + interference_sig = Fiber("fsimu",opt_signal.fs,"fiber_length",oneway_interference_meter*2/1000,"alpha",0.3,"D",0,"lambda0",1320,"gamma",0,"Dslope",0.07).process(opt_signal); + + interference_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",-sir).process(interference_sig); + +%%%%% Delay the "main" signal as in reality the interference is "older" than the main signal %%%%%% + [main_sig,dly] = opt_signal.delay("delay_meter",oneway_interference_meter*2); + + main_sig.power; + interference_sig.power; + +%%%%% ADD Interference and Main Signal %%%%%% + combined_sig = main_sig + interference_sig; + +%%%%% Cut (due to the delays there is a jump in the signals) %%%%%% + if dly == 0;dly = 1;end + combined_sig.signal = combined_sig.signal(ceil(dly):end); + +%%%%% Propagate through fiber %%%%%% + opt_signal = Fiber("fsimu",combined_sig.fs,"fiber_length",link_total_meter/1000,"alpha",0.3,"D",0,"lambda0",1320,"gamma",0,"Dslope",0.07).process(combined_sig); + +end \ No newline at end of file diff --git a/projects/MPI_August/mpi_uncorrelated.m b/projects/MPI_August/mpi_uncorrelated.m new file mode 100644 index 0000000..716fe91 --- /dev/null +++ b/projects/MPI_August/mpi_uncorrelated.m @@ -0,0 +1,206 @@ + +%% Parameter to simulate and save +params = struct; + +params.M = [4]; +params.datarate = [224]; +params.sir = [24]; %decibel = attenuation of interference path +params.laser_linewidth = [1e6]; +params.pn_key = [11]; +params.vbias_rel = [0.5]; +params.rop = -5; + +params.clipfactor = [1.5]; +params.rrcalpha = [0.1]; + +name = ['wh_',strrep(num2str(now),'.','')]; + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("rop"); +wh.addStorage("txpapr"); +wh.addStorage("er"); + +%% Init Params +link_length = 10000; %meter + +endcnt = prod(wh.dim); +cnt=0; + +disp(['Start Simulation of ',num2str(endcnt),' loops...']) +tic +for M = wh.parameter.M.values + for datarate = wh.parameter.datarate.values + for rrcalpha = wh.parameter.rrcalpha.values + + %% SETUP HERE: %% + kover = 8; + M8199 = M8199A("kover",kover); + fdac = M8199.fdac; + fsym = round(datarate / log2(M))*1e9; + + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + % MAIN SIGNAL + %%%%% Symbol Generation %%%%%% + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",2).process(); + %%%%% AWG %%%%%% + El_sig = M8199.process(Digi_sig); + El_sig = El_sig.*0.7222; + El_sig.signal = awgn(El_sig.signal,20,'measured',1); + %%%%% Lowpass before Modulator %%%%%% + El_sig = Filter('filtdegree',2,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",15).process(El_sig); + + + + % INTERFERENCE SIGNAL + %%%%% Symbol Generation %%%%%% + [Digi_sig_i,Symbols_i,Bits_i] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",1).process(); + %%%%% AWG %%%%%% + El_sig_i = M8199.process(Digi_sig_i); + El_sig_i = El_sig_i.*0.7222; + El_sig_i.signal = awgn(El_sig_i.signal,20,'measured',2); + %%%%% Lowpass before Modulator %%%%%% + El_sig_i = Filter('filtdegree',2,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_i); + El_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",15).process(El_sig_i); + + + for laser_linewidth = wh.parameter.laser_linewidth.values + for pn_key = wh.parameter.pn_key.values + for vbias_rel = wh.parameter.vbias_rel.values + + + + % MAIN SIGNAL + %%%%% MODULATE E/O CONVERSION %%%%%% + u_pi = 2.9; + vbias = -vbias_rel*u_pi; + [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key).process(El_sig); + Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + Opt_sig = Optfilter.process(Opt_sig); + + + % INTERFERENCE SIGNAL + %%%%% MODULATE E/O CONVERSION %%%%%% + u_pi = 2.9; + vbias = -vbias_rel*u_pi; + [Opt_sig_i] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key+1).process(El_sig_i); + Opt_sig_i = Optfilter.process(Opt_sig_i); + + + + j_ = wh.parameter.sir.length; + i_ = wh.parameter.rop.length; + ber=zeros(j_,i_); + patten=zeros(j_,i_); + + for j = 1:j_ + sir = wh.parameter.sir.values(j); + + + %%%%% Interference Signal Fiber Prop %%%%%% + Opt_sig_i = Fiber("fsimu",Opt_sig_i.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig_i); + + Opt_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",Opt_sig.power-sir).process(Opt_sig_i); + + %%%%% ADD Interference and Main Signal %%%%%% + Opt_sig = Opt_sig_i + Opt_sig; + + %%%%% Interference Signal Fiber Prop %%%%%% + Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + + + % % MPI Channel + % Opt = channel_model_mpi(Opt_sig,link_length,mpi_path,sir); + + % Receiver ROP curve + for i = 1:i_ + rop=wh.parameter.rop.values(i); + + % Set ROP + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); + patten(j,i) = Rx_sig.power; + + %%%%%% Square Law %%%%%% + Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20).process(Rx_sig); + + %%%%%% Lowpass PhDiode %%%%%% + Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig); + % Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201); + + + %%%%%% Scope %%%%%% + fadc = 256e9; + Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); + Scpe_sig.plot("displayname","SIgnal after Scope","fignum",999); + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + + %%%%% EQUALIZE %%%%%% + % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + % Eq = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",25,"dfe_order",2,"sps",2,"decide",1); + + % Eq = EQ("Ne",[25,3,3],"Nb",[0,0,0],"training_length",4096,"training_loops",4,"dd_loops",4,"K",2,"DCmu",0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",1); + + Eq = VNLE("epochs_tr",7,"epochs_dd",7,"len_tr",4096,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[25,3,3],"sps",2,"decide",0); + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + % EQ_sig.normalize("mode","rms").plot('fignum',23,'displayname','before eq') + + + %%%%% DEMAP %%%%%% + Rx_bits = PAMmapper(M,0).demap(EQ_sig); + + %%%% Look at Pam levels %%%%% + if 1 + a = PAMmapper(M,0).separate_pamlevels(EQ_sig); + figure(14);hold on;scatter(1:EQ_sig.length,a,1,'.'); + end + + % BER + [~,errors_bm,ber(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + cnt = cnt+1; + disp(['BER: ',sprintf('%.1E',ber(j,i)),' - - ROP: ',num2str(Rx_sig.power),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + + + end + end + + for j = 1:j_ + sir = wh.parameter.sir.values(j); + for i = 1:i_ + rop=wh.parameter.rop.values(i); + wh.addValueToStorage(ber(j,i),'ber',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); + wh.addValueToStorage(patten(j,i),'rop',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); + wh.addValueToStorage(er,'er',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); + end + + end + + + end + toc + disp(['Simulated: ',num2str(cnt/endcnt*100),' %']); + save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\',name,'.mat'],"wh"); + end + end + end + end +end + + + + + diff --git a/projects/MPI_August/rx_model.m b/projects/MPI_August/rx_model.m new file mode 100644 index 0000000..cef6a9e --- /dev/null +++ b/projects/MPI_August/rx_model.m @@ -0,0 +1,45 @@ +function rx_bits = rx_model(Rx,Ref,Phdiod,Lp_phdiod,Scpe,Eq,Pmap) + +%%%%% Local Parameter %%%%%%%% +fsym = Ref.fs; + + +%%%%%% Square Law %%%%%% + Rx_sig = Phdiod.process(Rx); + +%%%%%% Lowpass PhDiode %%%%%% + Rx_sig = Lp_phdiod.process(Rx_sig); + % Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201); + + +%%%%%% Scope %%%%%% + Scpe_sig = Scpe.process(Rx_sig); + Scpe_sig.plot("displayname","SIgnal after Scope","fignum",1999); + % Scpe_sig.spectrum("displayname",'after scope','fignum',123); + + +%%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",Scpe.fadc,"fs_out",2*fsym); + % Scpe_sig.plot('fignum',12345,'displayname','bla') + + +%%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Ref,"fs_ref",fsym); + + +%%%%% EQUALIZE %%%%%% + [EQ_sig] = Eq.process(Scpe_sig,Ref); + % EQ_sig.normalize("mode","rms").plot('fignum',23,'displayname','before eq') + + +%%%%% DEMAP %%%%%% + rx_bits = Pmap.demap(EQ_sig); + +%%%% Look at Pam levels %%%%% + if 1 + a = Pmap.separate_pamlevels(EQ_sig); + figure(14);hold on;scatter(1:EQ_sig.length,a,1,'.'); + end + +end + diff --git a/projects/MPI_August/tx_model.m b/projects/MPI_August/tx_model.m new file mode 100644 index 0000000..4d4e2c6 --- /dev/null +++ b/projects/MPI_August/tx_model.m @@ -0,0 +1,66 @@ + +function [el_sig,bits,symbol_seq] = tx_model(fsym,Pmap,Pform,Awg,options) + +arguments + fsym + Pmap + Pform + Awg + options.clipfactor +end + +%%%%% Local Parameter %%%%%%%% + + +%%%%% PRBS Generation in correct shape for Modulation Format %%%%%% + O = 17; %order of prbs + N = 2^(O-1); %length of prbs + [~,seed] = prbs(O,1); %initialize first seed of prbs + bitpattern=[]; + for i = 1:log2(Pmap.M) + [bitpattern(:,i),seed] = prbs(O,N,seed); + end + + if Pmap.M == 6 + bitpattern = reshape(bitpattern,[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); + end + + bits = Informationsignal(bitpattern); + +%%%%% Map to PAM %%%%%% + symbol_seq = Pmap.map(bits); + symbol_seq.fs = fsym; + + + +%%%%% Pulseforming %%%%%% +if 1 + X = Pform.process(symbol_seq); +else + X = symbol_seq; +end + +%%%%% Resample to f DAC %%%%%% + X = X.resample("fs_in",X.fs,"fs_out",Awg.fdac,"n",10,"beta",5); + +%%%%% Clip to PAM range %%%%%% +if 1 + min_ = min(symbol_seq.signal) * options.clipfactor ; + max_ = max(symbol_seq.signal) * options.clipfactor ; + X.signal = clip(X.signal,min_,max_); +end + +%%%%%Info about AWG input %%%%%% +pwr_ = X.power; % in dbm into 50 ohm +papr_ = papr(X.signal); +vpk_ = max(abs(X.signal)); +vswing_ = max(X.signal)-min(X.signal); +awg_string = ['Digital AWG Input: Pout: ',num2str(round(pwr_,2)),' dBm; PAPR: ',num2str(round(papr_,2)),'; Vpk: ',num2str(round(vpk_,2)),'; Vswing: ',num2str(round(vswing_,2)),'']; +% isp(awg_string); + +%%%%% AWG %%%%%% +el_sig = Awg.process(X); +el_sig= el_sig.*0.7222; + +end diff --git a/projects/MPI_Juni/EQ_optimization/eq_optimization.m b/projects/MPI_Juni/EQ_optimization/eq_optimization.m new file mode 100644 index 0000000..31060b4 --- /dev/null +++ b/projects/MPI_Juni/EQ_optimization/eq_optimization.m @@ -0,0 +1,99 @@ + +clear +Ref = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\EQ_optimization\ref_sig_with_mpi.mat",'Ref'); +Rx_sig = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\EQ_optimization\rx_sig_with_mpi.mat",'Rx_sig'); +Bits = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\EQ_optimization\ref_bits_prms.mat","Bits"); + +M = 4; +Ref = Ref.Ref; +Rx_sig = Rx_sig.Rx_sig; +Bits = Bits.Bits; + +kover = 8; +fdac = 256e9; +fadc = 160e9; + +Lp_phd = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + +Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + +Scp = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe); + +fadc = Scp.fadc; + +% SETUP TX Model +Pmap = PAMmapper(M,0); + +% SETUP RX Model +Phd = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); + +% Eq = EQ_silas_ofc("Ne",[50,2,2],"Nb",[2,0,0],"trainlength",4096,... +% "sps",2,... +% "mu_dc_dd",[0.5, 0.5, 0.5, 0.5],... +% "mu_dc_train",0.0,... +% "mu_ffe_train",0.00,... +% "mu_dfe_train",0.005,... +% "mu_ffe_dd",[0.0004 0.0004 0.0004],... +% "mu_dfe_dd",0.005,... +% "ddloops",3,... +% "trainloops",4,... +% "eq_parallelization_blocklength",0, ... +% "eq_updatelatency",1,... +% "eq_avg_blocklength",0); + + +if 0 + Eq = EQ("Ne",[21,4,4],"Nb",[0,0,0],"training_length",4096,"training_loops",4,"dd_loops",4,"K",2,"DCmu",0,"DDmu",[0.0004 0.0005 0.0006 0.0003 ],"DFEmu",0.000,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + + + % % RX Model + Rx_bits = rx_model(Rx_sig,Ref,Phd,Lp_phd,Scp,Eq,Pmap); + + % % BER + [~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + disp(['BER: ',sprintf('%.1E',ber),' - - ROP: ',num2str(Rx_sig.power),'dBm - - PAM-',num2str(M),' - - ']); +end + +%%%%% Local Parameter %%%%%%%% + fsym = Ref.fs; + +%%%%%% Square Law %%%%%% + Rx_sig = Phd.process(Rx_sig); + +%%%%%% Lowpass PhDiode %%%%%% + Rx_sig = Lp_phd.process(Rx_sig); + % Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201); + +%%%%%% Scope %%%%%% + Scpe_sig = Scp.process(Rx_sig); + % Scpe_sig.spectrum("displayname",'after scope','fignum',123); + +%%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",Scp.fadc,"fs_out",2*fsym); + % Scpe_sig.plot('fignum',12345,'displayname','bla') + +%%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Ref,"fs_ref",fsym); + + + +Scpe_sig = Scpe_sig.normalize("mode","rms"); + + +Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",20,"sps",2,"decide",1); +% Eq = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",21,"dfe_order",0,"sps",2,"decide",1); +% Eq = VNLE("epochs_tr",7,"epochs_dd",7,"len_tr",4096,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[21,5,5],"sps",2,"decide",1); + +Eq_sig = Eq.process(Scpe_sig,Ref); + +Rx_bits = Pmap.demap(Eq_sig); + +% BER +[~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + +disp(['BER: ',sprintf('%.1E',ber),' - - ROP: ',num2str(Rx_sig.power),'dBm - - PAM-',num2str(M),' - - ']); + diff --git a/projects/MPI_Juni/EQ_optimization/ref_bits_prms.mat b/projects/MPI_Juni/EQ_optimization/ref_bits_prms.mat new file mode 100644 index 0000000..4ebf4b8 Binary files /dev/null and b/projects/MPI_Juni/EQ_optimization/ref_bits_prms.mat differ diff --git a/projects/MPI_Juni/EQ_optimization/ref_sig_with_mpi.mat b/projects/MPI_Juni/EQ_optimization/ref_sig_with_mpi.mat new file mode 100644 index 0000000..7881c60 Binary files /dev/null and b/projects/MPI_Juni/EQ_optimization/ref_sig_with_mpi.mat differ diff --git a/projects/MPI_Juni/EQ_optimization/rx_sig_with_mpi.mat b/projects/MPI_Juni/EQ_optimization/rx_sig_with_mpi.mat new file mode 100644 index 0000000..a4ded7d Binary files /dev/null and b/projects/MPI_Juni/EQ_optimization/rx_sig_with_mpi.mat differ diff --git a/projects/MPI_Juni/ber_curve.m b/projects/MPI_Juni/ber_curve.m new file mode 100644 index 0000000..50d0815 --- /dev/null +++ b/projects/MPI_Juni/ber_curve.m @@ -0,0 +1,46 @@ +function ber_curve(wh,options) + arguments + wh + options.DisplayName = ''; + end + + M = wh.parameter.M.values(1); + datarate = wh.parameter.datarate.values; + sir = wh.parameter.sir.values(1); + laser_linewidth = wh.parameter.laser_linewidth.values(1); + pn_key = wh.parameter.pn_key.values(1); + vbias_rel = wh.parameter.vbias_rel.values(1); + rop = wh.parameter.rop.values; + clipfactor = wh.parameter.clipfactor.values; + rrcalpha = wh.parameter.rrcalpha.values(1); + + cols = linspecer(numel(datarate)); + + cnt = 0; + for dr = datarate + cnt = cnt+1; + bers = wh.getStoValue('ber',M,dr,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor, rrcalpha); + rops = wh.getStoValue('rop',M,dr,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor, rrcalpha); + ers = unique(wh.getStoValue('er',M,dr,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor, rrcalpha)); + + if 1 %isempty(options.DisplayName) + dn = ['M: ',num2str(M),'; Rate: ',num2str(dr),'; Vbias rel: ',num2str(vbias_rel), '; ER: ', num2str(ers)]; + else + dn = options.DisplayName; + end + + % Create the initial plot + figure(43); + hold on; % Retain the plot so new points can be added without complete redraw + plot(rops,bers,"LineWidth",1,"LineStyle","--","Marker",".","MarkerSize",10,"DisplayName",string(dn),'Color',cols(cnt,:)); + + end + yline(3.8e-3,'DisplayName','HD-FEC'); + xlabel('Received Optical Power (dBm)'); + ylabel('Bit Error Rate (BER)'); + title('Bit Error Rate vs. Received Optical Power'); + set(gca,'yscale','log'); + grid on; + legend + +end diff --git a/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_with_MPI.mat b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_with_MPI.mat new file mode 100644 index 0000000..b42a673 Binary files /dev/null and b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_with_MPI.mat differ diff --git a/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_with_MPI_linearEQ.mat b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_with_MPI_linearEQ.mat new file mode 100644 index 0000000..af84c79 Binary files /dev/null and b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_with_MPI_linearEQ.mat differ diff --git a/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_wo_MPI.mat b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_wo_MPI.mat new file mode 100644 index 0000000..61469f2 Binary files /dev/null and b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_wo_MPI.mat differ diff --git a/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_wo_MPI_linearEQ.mat b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_wo_MPI_linearEQ.mat new file mode 100644 index 0000000..e1a110f Binary files /dev/null and b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_clipping_bias_wo_MPI_linearEQ.mat differ diff --git a/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_withMPI.mat b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_withMPI.mat new file mode 100644 index 0000000..0313b02 Binary files /dev/null and b/projects/MPI_Juni/ber_vs_er/ber_vs_er_pam4_withMPI.mat differ diff --git a/projects/MPI_Juni/ber_vs_er/tx_signal_curve.m b/projects/MPI_Juni/ber_vs_er/tx_signal_curve.m new file mode 100644 index 0000000..a356372 --- /dev/null +++ b/projects/MPI_Juni/ber_vs_er/tx_signal_curve.m @@ -0,0 +1,44 @@ +function tx_signal_curve(wh,options) +arguments + wh + options.DisplayName = ''; +end + +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sirs = wh.parameter.sir.values(1); +mpi_len = wh.parameter.mpi_pathlen.values(1); +laser_linewidths = wh.parameter.laser_linewidth.values(1); +pn_key = wh.parameter.pn_key.values; +vbias_rel = wh.parameter.vbias_rel.values(1); +rop = wh.parameter.rop.values(1); +clipfactor = wh.parameter.clipfactor.values; +rrcalpha = wh.parameter.rrcalpha.values(1); +% Create the initial plot +figure(4); +cols = linspecer(6); +cnt = 0; + + +bers(:) = wh.getStoValue('ber',M,datarate,sirs,mpi_len,laser_linewidths,pn_key,vbias_rel,rop,clipfactor, rrcalpha); +ers(:) = wh.getStoValue('er',M,datarate,sirs,mpi_len,laser_linewidths,pn_key,vbias_rel,rop,clipfactor, rrcalpha); + +dn = ['M: ',num2str(M),'; Rate: ',num2str(datarate),'; SIR ',num2str(sirs), 'dB; lw: ',num2str(laser_linewidths.*1e-6),' MHz']; + +hold on; % Retain the plot so new points can be added without complete redraw +% scatter(clipfactor,bers,"LineWidth",1,"Marker",".","DisplayName",string(dn),"MarkerEdgeColor",cols(cnt,:),'HandleVisibility','off'); +yyaxis left +plot(clipfactor,bers,"LineWidth",1,"Marker",".","DisplayName",'BER'); +yline(3.8e-3,'DisplayName','HD-FEC','HandleVisibility','off'); +set(gca,'yscale','log'); +grid on; +ylim([1e-5,4e-2]) + +yyaxis right +plot(clipfactor,ers,"LineWidth",1,"Marker",".","DisplayName",'Extinction Ratio in dB'); +xlabel('Clip Factor'); +ylabel('Bit Error Rate (BER)'); +title(['Bit Error Rate vs. Clip Factor: ',dn]); + +legend +end diff --git a/projects/MPI_Juni/clip_optimization/pam4_clipping.fig b/projects/MPI_Juni/clip_optimization/pam4_clipping.fig new file mode 100644 index 0000000..647259a Binary files /dev/null and b/projects/MPI_Juni/clip_optimization/pam4_clipping.fig differ diff --git a/projects/MPI_Juni/clip_optimization/pam6_clipping.fig b/projects/MPI_Juni/clip_optimization/pam6_clipping.fig new file mode 100644 index 0000000..b345887 Binary files /dev/null and b/projects/MPI_Juni/clip_optimization/pam6_clipping.fig differ diff --git a/projects/MPI_Juni/clip_optimization/pam8_clipping.fig b/projects/MPI_Juni/clip_optimization/pam8_clipping.fig new file mode 100644 index 0000000..e98859f Binary files /dev/null and b/projects/MPI_Juni/clip_optimization/pam8_clipping.fig differ diff --git a/projects/MPI_Juni/clip_optimization/pam_4_6_8.fig b/projects/MPI_Juni/clip_optimization/pam_4_6_8.fig new file mode 100644 index 0000000..93227cb Binary files /dev/null and b/projects/MPI_Juni/clip_optimization/pam_4_6_8.fig differ diff --git a/projects/MPI_Juni/dataset_2506.mat b/projects/MPI_Juni/dataset_2506.mat new file mode 100644 index 0000000..a37e13e Binary files /dev/null and b/projects/MPI_Juni/dataset_2506.mat differ diff --git a/projects/MPI_Juni/extrcurve.m b/projects/MPI_Juni/extrcurve.m new file mode 100644 index 0000000..b377e57 --- /dev/null +++ b/projects/MPI_Juni/extrcurve.m @@ -0,0 +1,40 @@ +function ercurve(wh) + + M = wh.parameter.M.values(1); + datarate = wh.parameter.datarate.values(1); + sir = wh.parameter.sir.values(1); + laser_linewidth = wh.parameter.laser_linewidth.values(1); + pn_key = wh.parameter.pn_key.values(1); + vbias_rel = wh.parameter.vbias_rel.values; + rop = wh.parameter.rop.values; + clipfactor = wh.parameter.clipfactor.values; + rrcalpha = wh.parameter.rrcalpha.values(1); + + + for vbr = vbias_rel + + bers = wh.getStoValue('ber',M,datarate,sir,laser_linewidth,pn_key,vbr,rop,clipfactor, rrcalpha); + rops = wh.getStoValue('rop',M,datarate,sir,laser_linewidth,pn_key,vbr,rop,clipfactor, rrcalpha); + ers = unique(wh.getStoValue('er',M,datarate,sir,laser_linewidth,pn_key,vbr,rop,clipfactor, rrcalpha)); + + if 1 %isempty(options.DisplayName) + dn = ['M: ',num2str(M),'; Rate: ',num2str(datarate),'; Vbias rel: ',num2str(vbr), '; ER: ', num2str(ers)]; + else + dn = options.DisplayName; + end + + % Create the initial plot + figure(43); + hold on; % Retain the plot so new points can be added without complete redraw + plot(rops,bers,"LineWidth",1,"LineStyle","--","Marker",".","MarkerSize",10,"DisplayName",string(dn)); + + end + yline(3.8e-3,'DisplayName','HD-FEC'); + xlabel('Received Optical Power (dBm)'); + ylabel('Bit Error Rate (BER)'); + title('Bit Error Rate vs. Received Optical Power'); + set(gca,'yscale','log'); + grid on; + legend + +end \ No newline at end of file diff --git a/projects/MPI_Juni/linewidth_curve.m b/projects/MPI_Juni/linewidth_curve.m new file mode 100644 index 0000000..6d75000 --- /dev/null +++ b/projects/MPI_Juni/linewidth_curve.m @@ -0,0 +1,47 @@ +function linewidth_curve(wh,options) + arguments + wh + options.DisplayName = ''; + end + + M = wh.parameter.M.values(1); + datarate = wh.parameter.datarate.values(1); + sir = wh.parameter.sir.values(6); + mpi_len = wh.parameter.mpi_pathlen.values; + laser_linewidth = wh.parameter.laser_linewidth.values; + pn_key = wh.parameter.pn_key.values; + vbias_rel = wh.parameter.vbias_rel.values(1); + rop = wh.parameter.rop.values(1); + clipfactor = wh.parameter.clipfactor.values(1); + rrcalpha = wh.parameter.rrcalpha.values(1); + + cols = linspecer(numel(mpi_len)); + % Create the initial plot + figure(44); + hold on; % Retain the plot so new points can be added without complete redraw + for l = 1:numel(mpi_len) + + cnt = 0; + for i = 1:numel(laser_linewidth) + + bers(:,i) = wh.getStoValue('ber',M,datarate,sir,mpi_len(l),laser_linewidth(i),pn_key,vbias_rel,rop,clipfactor, rrcalpha); + %ers(:,i) = unique(wh.getStoValue('er',M,datarate,sir,mpi_len,laser_linewidth(i),pn_key,vbias_rel,rop,clipfactor, rrcalpha)); + + dn = ['Interference Path Length: ',num2str(mpi_len(l)),'m; Rate: ',num2str(datarate)]; + + end + + scatter(laser_linewidth.*1e-6,bers,"LineWidth",1,"Marker",".","DisplayName",string(dn),"MarkerEdgeColor",cols(l,:),'HandleVisibility','off'); + plot(laser_linewidth.*1e-6,mean(bers),"LineWidth",1,"Marker",".","DisplayName",string(dn),"MarkerEdgeColor",cols(l,:),"Color",cols(l,:)); + + end + ylim([1e-4, 1e-2]); + yline(3.8e-3,'DisplayName','HD-FEC'); + xlabel('Laser Linewidth in MHz'); + ylabel('Bit Error Rate (BER)'); + title(['Bit Error Rate vs. Interference Path Length @ SIR ',num2str(sir), 'dB']); + set(gca,'yscale','log'); + grid on; + legend + +end diff --git a/projects/MPI_Juni/main_simulation.m b/projects/MPI_Juni/main_simulation.m new file mode 100644 index 0000000..f8d46b2 --- /dev/null +++ b/projects/MPI_Juni/main_simulation.m @@ -0,0 +1,245 @@ + +%% Parameter to simulate and save +params = struct; + +params.M = [4]; +params.datarate = [224]; +params.sir = [25]; %decibel = attenuation of interference path +params.mpi_pathlen = [50]; +params.laser_linewidth = [1e6]; +params.pn_key = [15]; +params.vbias_rel = [0.5]; +params.rop = 0; + +params.clipfactor = [1.3]; +params.rrcalpha = [0.1]; + +name = ['wh_',strrep(num2str(now),'.','')]; + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("rop"); +wh.addStorage("txpapr"); +wh.addStorage("er"); + +%% Init Params +link_length = 10000; %meter + +endcnt = prod(wh.dim); +cnt=0; + +disp(['Start Simulation of ',num2str(endcnt),' loops...']) +tic +for M = wh.parameter.M.values + for datarate = wh.parameter.datarate.values + for rrcalpha = wh.parameter.rrcalpha.values + + %% SETUP HERE: %% + kover = 8; + Awg = M8199A("kover",kover); + fdac = Awg.fdac; + fsym = round(datarate / log2(M))*1e9; + + Lp_awg = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true); + Lp_mod = Filter('filtdegree',2,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true); + Lp_opt = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + Lp_phd = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + + fadc = 256e9; + Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + + Scp = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe); + + fadc = Scp.fadc; + + % figure(222) + % Lp_awg.showHere; + % Lp_mod.showHere; + % Lp_opt.showHere; + % Lp_phd.showHere; + % Lp_scpe.showHere; + + % SETUP TX Model + Pmap = PAMmapper(M,0); + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + % SETUP RX Model + Phd = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); + + % Eq = EQ_silas("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + % "sps",2,... + % "mu_dc_dd",0.1,... + % "mu_dc_train",0.0,... + % "mu_ffe_train",0.00,... + % "mu_dfe_train",0.005,... + % "mu_ffe_dd",[0.0004 0.0004 0.0004],... + % "mu_dfe_dd",0.005,... + % "ddloops",3,... + % "trainloops",4,... + % "eq_parallelization_blocklength",0, ... + % "eq_updatelatency",0,... + % "eq_avg_blocklength",1000); + + % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + + % Eq = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",25,"dfe_order",2,"sps",2,"decide",1); + + % Eq = EQ_silas_sliding_window_dc_removal("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + % "sps",2,... + % "mu_dc_dd",0.05,... + % "mu_dc_train",0.05,... + % "mu_ffe_train",0,... + % "mu_dfe_train",0.005,... + % "mu_ffe_dd",[0.0004 0.0004 0.0004],... + % "mu_dfe_dd",0.0004,... + % "ddloops",4,... + % "trainloops",4,... + % "eq_blocklength",0,... + % "eq_updatelatency",0); + + % Eq = EQ("Ne",[25,3,3],"Nb",[0,0,0],"training_length",4096,"training_loops",4,"dd_loops",4,"K",2,"DCmu",0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",1); + + Eq = VNLE("epochs_tr",7,"epochs_dd",7,"len_tr",4096,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[25,3,3],"sps",2,"decide",1); + + for clipfactor = wh.parameter.clipfactor.values + % TX model + [El,Bits,Ref] = tx_model(fsym,Pmap,Pform,Awg,"clipfactor",clipfactor); + + % Laser and Modulation + El = Lp_mod.process(El); + El = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",15).process(El); + pwr_ = El.power; % in dbm into 50 ohm + papr_ = papr(El.signal); + vpk_ = max(abs(El.signal)); + vswing_ = max(El.signal)-min(El.signal); + awg_string = ['AWG: Pout: ',num2str(round(pwr_,2)),' dBm (50 Ohm); PAPR: ',num2str(round(papr_,2)),'; Vpk: ',num2str(round(vpk_,2)),' V; Vswing: ',num2str(round(vswing_,2)),' V']; + disp(awg_string); + + for laser_linewidth = wh.parameter.laser_linewidth.values + for pn_key = wh.parameter.pn_key.values + for vbias_rel = wh.parameter.vbias_rel.values + + u_pi = 2.9; + vbias = -vbias_rel*u_pi; + + extmodlaser = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key); + % El = El.normalize("mode","oneone").*u_pi.*0.4; + [modulated,extmodlaser] = extmodlaser.process(El); + + if 1 + figure(10) + hold on + scatter(El.signal(1:100000)+vbias,(abs(modulated.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') + ylim([0 4]); + xlim([-u_pi/2, u_pi/2]+vbias); + xlabel('Input in V') + ylabel('abs(Output) in mW') + % Define properties + boxPosition = [0.15 0.86 0.2 0.05]; % Position for the first box [x y width height] + boxColor = [0.9 0.9 0.9]; % Light grey background color + boxEdgeColor = 'k'; % Black edge color + boxLineStyle = '--'; % Dashed line style + boxFontWeight = 'bold'; % Bold font + + % Create first annotation box for Power + annotation('textbox', boxPosition, ... + 'String', ['V pi: ',num2str(u_pi),' V'], ... + 'BackgroundColor', boxColor, ... + 'EdgeColor', boxEdgeColor, ... + 'LineStyle', boxLineStyle, ... + 'FontWeight', boxFontWeight, ... + 'HorizontalAlignment', 'center'); + + % Adjust position for the second box (slightly to the right) + boxPosition = [0.37 0.86 0.2 0.05]; % Adjusted position + + % Create second annotation box for PAPR + annotation('textbox', boxPosition, ... + 'String', ['V bias: ',num2str(vbias),' V'], ... + 'BackgroundColor', boxColor, ... + 'EdgeColor', boxEdgeColor, ... + 'LineStyle', boxLineStyle, ... + 'FontWeight', boxFontWeight, ... + 'HorizontalAlignment', 'center'); + + modulated.eye(fsym,M); + + modulated.spectrum("fignum",112,"displayname",'transmit spectrum'); + + modulated = Lp_opt.process(modulated); + + end + + er = modulated.extinctionratio(fsym,M); + + if 1 + + m_ = wh.parameter.mpi_pathlen.length; + j_ = wh.parameter.sir.length; + i_ = wh.parameter.rop.length; + ber=zeros(m_,j_,i_); + patten=zeros(m_,j_,i_); + + for m = 1:m_ + mpi_path = wh.parameter.mpi_pathlen.values(m); + for j = 1:j_ + sir = wh.parameter.sir.values(j); + + % MPI Channel + Opt = channel_model_mpi(modulated,link_length,mpi_path,sir); + + for i = 1:i_ + rop=wh.parameter.rop.values(i); + + % Set ROP + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt); + patten(m,j,i) = Rx_sig.power; + + % RX Model + Rx_bits = rx_model(Rx_sig,Ref,Phd,Lp_phd,Scp,Eq,Pmap); + + % BER + [~,errors_bm,ber(m,j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + cnt = cnt+1; + disp(['BER: ',sprintf('%.1E',ber(m,j,i)),' - - ROP: ',num2str(Rx_sig.power),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + + end + end + end + + for m = 1:m_ + mpi_path = wh.parameter.mpi_pathlen.values(m); + for j = 1:j_ + sir = wh.parameter.sir.values(j); + for i = 1:i_ + rop=wh.parameter.rop.values(i); + wh.addValueToStorage(ber(m,j,i),'ber',M,datarate,sir,mpi_path,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); + wh.addValueToStorage(patten(m,j,i),'rop',M,datarate,sir,mpi_path,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); + wh.addValueToStorage(er,'er',M,datarate,sir,mpi_path,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha); + end + end + end + + + end + + end + toc + disp(['Simulated: ',num2str(cnt/endcnt*100),' %']); + save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\',name,'.mat'],"wh"); + end + end + end + end + end +end + + + + + diff --git a/projects/MPI_Juni/mpilength_curve.m b/projects/MPI_Juni/mpilength_curve.m new file mode 100644 index 0000000..e13f5ce --- /dev/null +++ b/projects/MPI_Juni/mpilength_curve.m @@ -0,0 +1,56 @@ +function mpilength_curve(wh,options) +arguments + wh + options.DisplayName = ''; +end + +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(2); +sirs = wh.parameter.sir.values(1); +mpi_len = wh.parameter.mpi_pathlen.values; +laser_linewidths = wh.parameter.laser_linewidth.values; +pn_key = wh.parameter.pn_key.values; +vbias_rel = wh.parameter.vbias_rel.values(1); +rop = wh.parameter.rop.values(1); +clipfactor = wh.parameter.clipfactor.values(1); +rrcalpha = wh.parameter.rrcalpha.values(1); +% Create the initial plot +figure(4); +cols = linspecer(6); +cnt = 0; + +for lw = 1:numel(laser_linewidths) + laser_linewidth = laser_linewidths(lw); + n=ceil(sqrt(wh.parameter.laser_linewidth.length)); + %subplot(n,n,sp); + cnt=cnt+6; + for sir = sirs + + for i = 1:numel(mpi_len) + + bers(:,i) = wh.getStoValue('ber',M,datarate,sir,mpi_len(i),laser_linewidth,pn_key,vbias_rel,rop,clipfactor, rrcalpha); + % ers(:,i) = unique(wh.getStoValue('er',M,datarate,sir,mpi_len(i),laser_linewidth,pn_key,vbias_rel,rop,clipfactor, rrcalpha)); + + dn = ['M: ',num2str(M),'; Rate: ',num2str(datarate),'; SIR ',num2str(sir), 'dB; lw: ',num2str(laser_linewidth.*1e-6),' MHz']; + + end + + + hold on; % Retain the plot so new points can be added without complete redraw + scatter(mpi_len,bers,"LineWidth",1,"Marker",".","DisplayName",string(dn),"MarkerEdgeColor",cols(cnt,:),'HandleVisibility','off'); + plot(mpi_len,mean(bers),"LineWidth",1,"Marker",".","DisplayName",string(dn),"MarkerEdgeColor",cols(cnt,:),"Color",cols(cnt,:)); + + + end + + yline(3.8e-3,'DisplayName','HD-FEC'); + xlabel('MPI Path in Meter'); + ylabel('Bit Error Rate (BER)'); + title(['Bit Error Rate vs. Interference Path Length; Linewidth: ',num2str(laser_linewidth.*1e-6), ' MHz']); + set(gca,'yscale','log'); + grid on; + ylim([1e-5,4e-2]) + +end +legend +end diff --git a/projects/MPI_Juni/sir_curve.m b/projects/MPI_Juni/sir_curve.m new file mode 100644 index 0000000..c5d4278 --- /dev/null +++ b/projects/MPI_Juni/sir_curve.m @@ -0,0 +1,46 @@ +function sir_curve(wh,options) + arguments + wh + options.DisplayName = ''; + end + + M = wh.parameter.M.values(1); + datarate = wh.parameter.datarate.values(1); + sir = wh.parameter.sir.values; + laser_linewidth = wh.parameter.laser_linewidth.values(1); + pn_key = wh.parameter.pn_key.values; + vbias_rel = wh.parameter.vbias_rel.values(1); + rop = wh.parameter.rop.values(1); + clipfactor = wh.parameter.clipfactor.values; + rrcalpha = wh.parameter.rrcalpha.values(1); + + cols = linspecer(numel(pn_key)); + + cnt = 0; + for pnk = pn_key + + cnt = cnt+1; + bers = wh.getStoValue('ber',M,datarate,sir,laser_linewidth,pnk,vbias_rel,rop,clipfactor, rrcalpha); + ers = unique(wh.getStoValue('er',M,datarate,sir,laser_linewidth,pnk,vbias_rel,rop,clipfactor, rrcalpha)); + + if 1 %isempty(options.DisplayName) + dn = ['M: ',num2str(M),'; Rate: ',num2str(datarate),'; Vbias rel: ',num2str(vbias_rel), '; ER: ', num2str(ers)]; + else + dn = options.DisplayName; + end + + % Create the initial plot + figure(43); + hold on; % Retain the plot so new points can be added without complete redraw + plot(sir,bers',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName",string(dn),'Color',cols(cnt,:)); + + end + yline(3.8e-3,'DisplayName','HD-FEC'); + xlabel('Signal to Interference Ratio (dB)'); + ylabel('Bit Error Rate (BER)'); + title('Bit Error Rate vs. SIR (MPI)'); + set(gca,'yscale','log'); + grid on; + legend + +end diff --git a/projects/MPI_Juni/tx_signal_curve.m b/projects/MPI_Juni/tx_signal_curve.m new file mode 100644 index 0000000..2666c6c --- /dev/null +++ b/projects/MPI_Juni/tx_signal_curve.m @@ -0,0 +1,53 @@ +function tx_signal_curve(wh,options) +arguments + wh + options.DisplayName = ''; +end + +M = wh.parameter.M.values(1); +datarate = wh.parameter.datarate.values(1); +sirs = wh.parameter.sir.values(1); +mpi_len = wh.parameter.mpi_pathlen.values(1); +laser_linewidths = wh.parameter.laser_linewidth.values(1); +pn_key = wh.parameter.pn_key.values; +vbias_rel = wh.parameter.vbias_rel.values; +rop = wh.parameter.rop.values(1); +clipfactor = wh.parameter.clipfactor.values; +rrcalpha = wh.parameter.rrcalpha.values(1); + +% Create the initial plot +figure(4); +cols = linspecer(numel(vbias_rel)); +cnt = 0; + +for vb = vbias_rel + + cnt = cnt+1; + bers(:) = wh.getStoValue('ber',M,datarate,sirs,mpi_len,laser_linewidths,pn_key,vb,rop,clipfactor, rrcalpha); + ers(:) = wh.getStoValue('er',M,datarate,sirs,mpi_len,laser_linewidths,pn_key,vb,rop,clipfactor, rrcalpha); + dn = ['M: ',num2str(M),'; Rate: ',num2str(datarate),'; SIR ',num2str(sirs), 'dB; lw: ',num2str(laser_linewidths.*1e-6),' MHz']; + + hold on; % Retain the plot so new points can be added without complete redraw + % scatter(clipfactor,bers,"LineWidth",1,"Marker",".","DisplayName",string(dn),"MarkerEdgeColor",cols(cnt,:),'HandleVisibility','off'); + plot(clipfactor,bers,"LineWidth",1,"Marker","x",'LineStyle','--',"DisplayName",['Vbias factor: ',num2str(vb)],"Color",cols(cnt,:)); + yline(3.8e-3,'DisplayName','HD-FEC','HandleVisibility','off'); + set(gca,'yscale','log'); + grid on; + ylim([1e-5,5e-1]) + + xlabel('Clip Factor'); + ylabel('Bit Error Rate (BER)'); + title(['Bit Error Rate vs. Clip Factor: ',dn]); + + legend +end + + +% figure(10101) +% y = vbias_rel; +% x = clipfactor; +% [X,Y] = meshgrid(x,y); +% Z = bers; +% contourf(X,Y,Z,30,'LineStyle','none'); + +end diff --git a/projects/MPI_Juni/wh_first_run.mat b/projects/MPI_Juni/wh_first_run.mat new file mode 100644 index 0000000..63735cf Binary files /dev/null and b/projects/MPI_Juni/wh_first_run.mat differ diff --git a/projects/MPI_analysis_Mai24/base_simulation.m b/projects/MPI_analysis_Mai24/base_simulation.m new file mode 100644 index 0000000..e6579dc --- /dev/null +++ b/projects/MPI_analysis_Mai24/base_simulation.m @@ -0,0 +1,124 @@ + +M=4; + +fdac = 256e9;%fsym; +fadc = 256e9; +fsym = [96:16:256].*1e9; + +%fsym = 160e9; +% 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)); +end + +% 2 ) Build Inf. signal class +bits = Informationsignal(bitpattern); + +% 3) Digi modulation -> PAM-M signal +digimod_out = PAMmapper(M,0).map(bits); + +% 5) AWG (lowpass, quantization, sample and hold) +kover = 8; +LP_awg = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*kover,"filterType",filtertypes.butterworth); + +powerlist = []; +for i = length(fsym):-1:1 + digimod_out.fs = fsym(i); + + X = Pulseformer("fsym",fsym(i),"fdac",fdac,"pulse","rrc","pulselength",16,"rrcalpha",0.1).process(digimod_out); + %X = digimod_out; + + X = AWG("fdac",fdac,"dac_min",-1,"dac_max",1,"lpf_active",1,"H_lpf",LP_awg,"kover",kover,"bit_resolution",5.5,"normalize2dac",1,"upsampling_method","samplehold").process(X); + + + % 6) Lowpass behavior before laser + X = LP_modulator.process(X); + + % 7) Normalize signal + X = X.normalize("mode","oneone"); + + % 1) Laser; Modulation -> OPTICAL DOMAIN + u_pi = 2; + vbias = -vb(m); + extmodlaser = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",X.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth(l),"randomkey",pn_key(pnk)); + E = X.*vp(n); + + [Opt,extmodlaser] = extmodlaser.process(E); + + figure(m) + hold on + scatter(E.signal(1:100000),(abs(Opt.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') + xlabel('Input in V') + ylabel('abs(Output) in mW') + + % ER = 10*log10(max(abs(Opt.signal).^2)/min(abs(Opt.signal).^2)); + + Opt = LP_opt.process(Opt); + + cspr(s,l,pnk,n,m) = Opt.cspr; + mod_out_pow(s,l,pnk,n,m) = Opt.power; + + % 2) ping pong fiber propagation + Interference_sig = Fiber("fsimu",Opt.fs,"fiber_length",mpi_path*2/1000,"alpha",0,"D",0,"lambda0",1310,"gamma",0).process(Opt); + + Interference_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",-sir(s)).process(Interference_sig); + + % In the meantime: delay the main signal + [Main_sig,dly] = Opt.delay("delay_meter",mpi_path*2); + + % Add + Combined_sig = Main_sig + Interference_sig; + + % Cut (due to the delays there is a jump in the signals) + if dly == 0;dly = 1;end + Combined_sig.signal = Combined_sig.signal(ceil(dly):end); + + % 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); + + + + + powerlist(i)=X.power; + + % Sample to 2x fsym + X = X.resample("fs_in",kover*fdac,"fs_out",2*fsym(i)); + + % Sync Rx signal with reference + [X,D,cuts] = X.tsynch("reference",digimod_out,"fs_ref",fsym(i)); + + [EQ_sig,EQ_sym] = EQ_silas("Ne",[50,0,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",0.00,... + "mu_dc_train",0.0,... + "mu_ffe_train",0.00,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0006 0.0003],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",3,... + "eq_parallelization_blocklength",1, ... + "eq_updatelatency",1,... + "eq_avg_blocklength",0).process(X,digimod_out); + + % Demap + Rx_Bits = PAMmapper(M,0).demap(EQ_sig); + + % BER + [~,errors_bm,BER(i),errors] = calc_ber(Rx_Bits.signal,bitpattern,"skip_front",100,"skip_end",150,"returnErrorLocation",1); +end + + +figure() +stem(fsym.*1e-9,powerlist); +xticks(fsym.*1e-9) \ No newline at end of file diff --git a/projects/MPI_analysis_Mai24/eq_compare.m b/projects/MPI_analysis_Mai24/eq_compare.m new file mode 100644 index 0000000..857d681 --- /dev/null +++ b/projects/MPI_analysis_Mai24/eq_compare.m @@ -0,0 +1,247 @@ + +%clear; +col = linspecer(6); + +% GENERATE SIR CURVE +if ismac + foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km'; +else + foldername = 'C:\Users\Silas\Documents\MATLAB\Datensätze\Labor_Datensatz_PAM4_MPI_OFC2023\pam4_10km_1km'; +end + +allfiles = dir(foldername); + +eq_updatelatency = 2; +eq_avg_blocklength =0; +eq_parallelization_blocklength = 92; +mudc = 0.01; + + +eq_ofc = EQ_silas_ofc("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",mudc,... + "mu_dc_train",mudc,... + "mu_ffe_train",0,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0004 0.0004],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",4,... + "eq_parallelization_blocklength",eq_parallelization_blocklength, ... + "eq_updatelatency",eq_updatelatency,... + "eq_avg_blocklength",eq_avg_blocklength); + +eq_normal = EQ_silas_ofc("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",0,... + "mu_dc_train",0,... + "mu_ffe_train",0,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0004 0.0004],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",4,... + "eq_parallelization_blocklength",eq_parallelization_blocklength, ... + "eq_updatelatency",eq_updatelatency,... + "eq_avg_blocklength",eq_avg_blocklength); + +mudc = 0.01; +eq_ff = EQ_silas("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",mudc,... + "mu_dc_train",mudc,... + "mu_ffe_train",0,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0004 0.0004],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",4,... + "eq_parallelization_blocklength",eq_parallelization_blocklength, ... + "eq_updatelatency",eq_updatelatency,... + "eq_avg_blocklength",0); + +% [ber_nml,sir_nml,fsym_nml] = sweepSIR(allfiles,eq_normal); +% [ber_ofc,sir_ofc,fsym_ofc] = sweepSIR(allfiles,eq_ofc); +% [ber_ff,sir_ff,fsym_ff] = sweepSIR(allfiles,eq_ff); + +% plotBerCurve(ber_nml,sir_nml,fsym_nml); +% plotBerCurve(ber_ofc,sir_ofc,fsym_ofc); +% plotBerCurve(ber_ff,sir_ff,fsym_ff); + +measurementpath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\Labor_Datensatz_PAM4_MPI_OFC2023\pam4_10km_1km\pam4__loop_10_92Gbd_19092023_1508.mat'; +% measurementpath = findSirAndFsym(allfiles,-30,92e9); +measurementpath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\Labor_Datensatz_PAM4_MPI_OFC2023\pam4_10km_1km\pam4__loop_30_92Gbd_19092023_1547.mat'; + +% comparePSD(measurementpath) + +% [totalbits,errors,ber_,loc,sir,fsym] = singleBerRun(measurementpath, eq_normal); +% [totalbits,errors,ber_,loc,sir,fsym] = singleBerRun(measurementpath, eq_ofc); + + +[~,~,ber_,~,sir,fsym] = singleBerRun(measurementpath, eq_ff); + + +mudc_loop = [0, 0.0001,0.0005, 0.001,0.005, 0.01:0.01:0.1, 0.2:0.1:1]; + +blocklength = [20:20:200]; + +parfor m = 1:length(mudc_loop) + eq_ff = EQ_silas("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",mudc_loop(m),... + "mu_dc_train",mudc_loop(m),... + "mu_ffe_train",0,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0004 0.0004],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",4,... + "eq_parallelization_blocklength",eq_parallelization_blocklength, ... + "eq_updatelatency",eq_updatelatency,... + "eq_avg_blocklength",100);%blocklength(m)); + + [~,~,ber_(m),~,sir(m),fsym(m)] = singleBerRun(measurementpath, eq_ff); +end + +figure(2) +hold on +plot(mudc_loop,ber_,'Marker','o') +ylim([1e-4, 1e-2]) +yline(3.8e-3); +set(gca, 'YScale', 'log'); +set(gca, 'XScale', 'log'); + + + +function cur_path = findSirAndFsym(allfiles,sir_desired,fsym_desired) + for i = 1:length(allfiles) + + if allfiles(i).bytes ~= 0 + cur_path = [allfiles(i).folder,filesep, allfiles(i).name]; + recorded_data = load(cur_path); + delete(findobj('Type','figure','Name','prms_compare')); + else + continue + end + + sir = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + fsym = recorded_data.saveStructTemp.common.f_sym; + + if abs(sir-sir_desired)<1 && abs(fsym-fsym_desired)<1e9 + return + end + end +end + + +function comparePSD(filepath) + recorded_data = load(filepath); + delete(findobj('Type','figure','Name','prms_compare')); + sir = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + fsym = recorded_data.saveStructTemp.common.f_sym; + + spectrum_plot(recorded_data.saveStructTemp.dp_tsynch_out',2*fsym); + + x = recorded_data.saveStructTemp.dp_tsynch_out'; + + pwelch(x,hamming(length(x)),[],length(x),2*fsym,"centered","power"); + + periodogram(x,hamming(length(x)),length(x),"centered","psd"); + +end + +function [totalbits,errors,ber_,loc,sir,fsym] = singleBerRun(filepath, eq_object) + + recorded_data = load(filepath); + delete(findobj('Type','figure','Name','prms_compare')); + sir = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + fsym = recorded_data.saveStructTemp.common.f_sym; + + % 0) build RX Signal + y_rx = Electricalsignal(recorded_data.saveStructTemp.dp_tsynch_out'); + y_rx.fs = 2.*fsym; + + % 0) build tx reference Signal for eq training + y_digimod = Informationsignal(recorded_data.saveStructTemp.digi_mod_out'); + y_digimod.fs = fsym; + + bits_ref = recorded_data.saveStructTemp.prms_out; + + [totalbits,errors,ber_,loc] = runPostProc(eq_object,y_rx,y_digimod,bits_ref); + + disp([class(eq_object),': fsym: ',num2str(fsym*1e-9),'GBd, SIR: ',num2str(sir),' ->> BER: ',sprintf('%2E',ber_)]); + +end + +function plotBerCurve(ber_,sir,fsym) + + rate = [92e9, 56e9]; + figure(95) + + for r = 1:length(rate) + + sorted = sortrows([sir(fsym == rate(r)); ber_(fsym == rate(r))]',1)'; + hold on + plot(abs(sorted(1,:)),sorted(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4;'],'LineWidth',1,'Marker','o','MarkerSize',5,'LineStyle','-','HandleVisibility','on'); + + + end + set(gca, 'YScale', 'log'); + yline(3.8e-3,'LineWidth',1, 'LineStyle','--','HandleVisibility','off'); + xlim([15,35]); + xlabel('SIR in dB'); + ylabel('BER'); +end + +function [ber_,sir,fsym] = sweepSIR(allfiles,eq_object) + + parfor i = 1:length(allfiles) + + if allfiles(i).bytes ~= 0 + recorded_data = load([allfiles(i).folder,filesep, allfiles(i).name]); + delete(findobj('Type','figure','Name','prms_compare')); + else + continue + end + + sir(i) = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + fsym(i) = recorded_data.saveStructTemp.common.f_sym; + + fdac(i) = recorded_data.saveStructTemp.common.f_DAC; + fadc(i) = recorded_data.saveStructTemp.common.f_ADC; + + % 0) build RX Signal + y_rx = Electricalsignal(recorded_data.saveStructTemp.dp_tsynch_out'); + y_rx.fs = 2.*fsym(i); + + % 0) build tx reference Signal for eq training + y_digimod = Informationsignal(recorded_data.saveStructTemp.digi_mod_out'); + y_digimod.fs = fsym(i); + + bits_ref = recorded_data.saveStructTemp.prms_out; + + [totalbits,errors,ber_(i),loc] = runPostProc(eq_object,y_rx,y_digimod,bits_ref); + + disp(['fsym: ',num2str(fsym(i)*1e-9),'GBd, SIR: ',num2str(sir(i)),' ->> BER: ',sprintf('%2E',ber_(i))]); + + end + +end + +function [totalbits,errors,ber_,loc] = runPostProc(eq_object,y_rx,y_ref,bits_ref) + + % 1) normlaize + y_rx = y_rx.normalize("mode","rms"); + + [Eq_out] = eq_object.process(y_rx,y_ref); + + % 3) digital demodulation object + digimod = PAMmapper(4,0); + + %estimated/ equalized symbol sequence + d_estimated = digimod.demap(Eq_out); + + % 4) BER calculation + [totalbits,errors,ber_,loc] = calc_ber(d_estimated.signal(1:end,:),bits_ref(1:end,:)',"skip_front",100,"skip_end",150,"returnErrorLocation",1); + +end \ No newline at end of file diff --git a/projects/MPI_analysis_Mai24/ofc_sim.m b/projects/MPI_analysis_Mai24/ofc_sim.m new file mode 100644 index 0000000..3e80085 --- /dev/null +++ b/projects/MPI_analysis_Mai24/ofc_sim.m @@ -0,0 +1,462 @@ + +clear; +col = linspecer(6); + +% GENERATE SIR CURVE +if ismac + foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km'; +else + foldername = 'C:\Users\Silas\Documents\MATLAB\Datensätze\Labor_Datensatz_PAM4_MPI_OFC2023\pam4_10km_1km'; +end + + +% SETTINGS +optimize_mudc = 0; +run_sir_sweep = 1; +run_feed_forward = 0; +run_baseline = 0; +run_ideal_dc_tap = 0; +plot_timesignal = 0; + +block_loop = [92]; + +for bl = 1:numel(block_loop) + + eq_parallelization_blocklength = block_loop(bl); + eq_updatelatency = 3; + eq_avg_blocklength =0; + + % FIND BEST MU DC + + if optimize_mudc + mudc_loop = [0, 0.0001,0.0005, 0.001,0.005, 0.01:0.01:0.1, 0.2:0.1:1]; + + current_filename = 'pam4__loop_14_92Gbd_19092023_1513.mat'; + % current_filename = 'pam4__loop_29_56Gbd_19092023_1546.mat'; + recorded_data = load([foldername,filesep, current_filename]); + + [best_mudc,ber]= optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop); + + figure(16) + hold on + scatter(mudc_loop(mudc_loop~=best_mudc),ber(mudc_loop~=best_mudc),'Marker','*','LineWidth',1); + set(gca, 'YScale', 'log'); + set(gca, 'XScale', 'log'); + scatter(best_mudc,ber(mudc_loop==best_mudc),100,'Marker','*','LineWidth',2) + yline(ber(1)); + xlim([0,1]) + end + + if run_sir_sweep + + %mudc = 0.04;%best_mudc; + mudc = 0.005; + + [ber,sir,fsym] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc); + + plotSirSweep(ber,sir,fsym,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + crossing = calculateCrossing([92e9, 56e9],ber,sir,fsym); + + req_sir_56(bl) = crossing(2,2) ; + req_sir_92(bl) = crossing(2,1) ; + + end + + +end + +if run_feed_forward + + current_filename = 'pam4__loop_14_92Gbd_19092023_1513.mat'; + recorded_data = load([foldername,filesep, current_filename]); + eq_avg_blocklength = [50,100,1000,3000]; + [best_block,ber]= optimizeAvgBlocklength(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,0); + + + best_block = 400; + + [ber_base,sir_base,fsym_base] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,best_block,0); + + plotSirSweep(ber_base,sir_base,fsym_base,eq_parallelization_blocklength,eq_updatelatency,best_block,0) + + crossing = calculateCrossing([92e9, 56e9],ber_base,sir_base,fsym_base); + + [ber_base,sir_base,fsym_base] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,0,0); + + plotSirSweep(ber_base,sir_base,fsym_base,eq_parallelization_blocklength,eq_updatelatency,0,0) + + crossing = calculateCrossing([92e9, 56e9],ber_base,sir_base,fsym_base); + +end + + +if run_baseline + + % baseline means no DC removal alg! + mudc = 0; + [ber_base,sir_base,fsym_base] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc); + + plotSirSweep(ber_base,sir_base,fsym_base,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + crossing = calculateCrossing([92e9, 56e9],ber_base,sir_base,fsym_base); + + req_sir_56_baseline = crossing(2,2) ; + req_sir_92_baseline = crossing(2,1) ; + +end + +if run_ideal_dc_tap + + current_filename = 'pam4__loop_29_56Gbd_19092023_1546.mat'; + recorded_data = load([foldername,filesep, current_filename]); + [mudc,~]= optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop); + + eq_parallelization_blocklength = 1; + eq_updatelatency = 1; + [ber_ideal,sir_ideal,fsym_ideal] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc); + + plotSirSweep(ber_ideal,sir_ideal,fsym_ideal,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + crossing = calculateCrossing([92e9, 56e9],ber_ideal,sir_ideal,fsym_ideal); + + req_sir_56_ideal = crossing(2,2) ; + req_sir_92_ideal = crossing(2,1) ; + +end + + + +if plot_timesignal + + current_filename = 'pam4__loop_26_92Gbd_19092023_1542.mat'; + recorded_data = load([foldername,filesep, current_filename]); + + sir = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + + fsym = 92e9; + eq_parallelization_blocklength = 1; + eq_updatelatency = 1; + eq_avg_blocklength = 0; + + % TRUE SYMBOLS + correct_symbols = recorded_data.saveStructTemp.digi_mod_out'; + + %3) PLAIN RECEIVED SIGNAL + disp('RX Signal') + y_rx = Electricalsignal(recorded_data.saveStructTemp.dp_tsynch_out'); + y_rx.fs = 2.*fsym; + rx_symbols = y_rx.resample("fs_in",2.*fsym,"fs_out",fsym); + rx_symbols = rx_symbols.normalize("mode","rms").signal; + disp(std(rx_symbols)); + scatterleveldependent(rx_symbols,correct_symbols,fsym); + + %1) PLOT WITH WITH ALGORITHM + if 0 + mudc_loop = [0, 0.0001,0.0005, 0.001,0.005, 0.01:0.01:0.1, 0.2:0.1:1]; + [mudc,~]= optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop); + end + + disp('ALGORITHM ') + mudc = 0.05; + results_alg = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + disp(results_alg.ber); + rx_symbols = results_alg.EQ_out.signal; + disp(std(rx_symbols)); + scatterleveldependent(rx_symbols,correct_symbols,fsym); + + %2) JUST EQ + disp('Just EQ') + mudc = 0; + results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + disp(results.ber) + rx_symbols = results.EQ_out.signal; + disp(std(rx_symbols)); + scatterleveldependent(rx_symbols,correct_symbols,fsym); + + + + +end + + +function plotSirSweep(ber,sir,fsym,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + + rate = [92e9, 56e9]; + figure(95) + + for r = 1:length(rate) + + sorted = sortrows([sir(fsym == rate(r)); ber(fsym == rate(r))]',1)'; + hold on + plot(abs(sorted(1,:)),sorted(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; mudc: ',num2str(mudc)],'LineWidth',1,'Marker','o','MarkerSize',5,'LineStyle','-','HandleVisibility','on'); + + end + + set(gca, 'YScale', 'log'); + yline(3.8e-3,'LineWidth',1, 'LineStyle','--','HandleVisibility','off'); + xlim([15,35]); + xlabel('SIR in dB'); + ylabel('BER'); + +end + +function [crossing] = calculateCrossing(rate,ber,sir,fsym) + + for r = 1:length(rate) + sorted = sortrows([sir(fsym == rate(r)); ber(fsym == rate(r))]',1)'; + berVals = sorted(2,:); + xAxis = sorted(1,:); + + hdfec = 3.8e-3 .* ones(1,length(sorted)); + crossing(:,r) = InterX([hdfec(:)';xAxis],[berVals;xAxis]) ; + end + +end + +function [ber,sir,fsym] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) + + % GENERATE SIR CURVE + if ismac + foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km'; + else + foldername = 'C:\Users\Silas\Documents\MATLAB\Datensätze\Labor_Datensatz_PAM4_MPI_OFC2023\pam4_10km_1km'; + end + + allfiles = dir(foldername); + + for i = 1:length(allfiles) + + if allfiles(i).bytes ~= 0 + current_filename = allfiles(i).name; + recorded_data = load([foldername,filesep, current_filename]); + else + continue + end + + results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + ber(i) = results.ber; + fsym(i) = results.fsym; + sir(i) = results.sir; + + disp(['fsym: ',num2str(results.fsym*1e-9),'GBd, SIR: ',num2str(results.sir),' ->> BER: ',sprintf('%2E',results.ber)]); + + end +end + +function [best_blocklength,ber] = optimizeAvgBlocklength(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength_loop,mudc) + + for j = 1:length(eq_avg_blocklength_loop) + eq_avg_blocklength = eq_avg_blocklength_loop(j); + + + results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + ber(j) = results.ber; + + disp(['fsym: ',num2str(results.fsym*1e-9),'GBd, SIR: ',num2str(results.sir),' ->> BER: ',sprintf('%2E',results.ber)]); + + end + + [~,pos]=min(ber); + best_blocklength = eq_avg_blocklength_loop(pos); + +end + +function [best_mudc,ber] = optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop) + + parfor j = 1:length(mudc_loop) + mudc = mudc_loop(j); + + + results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); + ber(j) = results.ber; + + disp(['fsym: ',num2str(results.fsym*1e-9),'GBd, SIR: ',num2str(results.sir),' ->> BER: ',sprintf('%2E',results.ber)]); + + end + + [~,pos]=min(ber); + best_mudc = mudc_loop(pos); + +end + +function results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength) + + results.sir = recorded_data.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + results.fsym = recorded_data.saveStructTemp.common.f_sym; + + results.fdac = recorded_data.saveStructTemp.common.f_DAC; + results.fadc = recorded_data.saveStructTemp.common.f_ADC; + + + % 0) build RX Signal + y_rx = Electricalsignal(recorded_data.saveStructTemp.dp_tsynch_out'); + y_rx.fs = 2.*results.fsym; + + % 0) build tx reference Signal for eq training + y_digimod = Informationsignal(recorded_data.saveStructTemp.digi_mod_out'); + y_digimod.fs = results.fsym; + + + % 1) normlaize + y_rx = y_rx.normalize("mode","rms"); + + eq = EQ_silas_ofc("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + "sps",2,... + "mu_dc_dd",mudc,... + "mu_dc_train",mudc,... + "mu_ffe_train",0,... + "mu_dfe_train",0.005,... + "mu_ffe_dd",[0.0004 0.0004 0.0004],... + "mu_dfe_dd",0.005,... + "ddloops",3,... + "trainloops",4,... + "eq_parallelization_blocklength",eq_parallelization_blocklength, ... + "eq_updatelatency",eq_updatelatency,... + "eq_avg_blocklength",eq_avg_blocklength); + % + % eq = EQ_silas("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + % "sps",2,... + % "mu_dc_dd",mudc,... + % "mu_dc_train",mudc,... + % "mu_ffe_train",0,... + % "mu_dfe_train",0.005,... + % "mu_ffe_dd",[0.0004 0.0004 0.0004],... + % "mu_dfe_dd",0.005,... + % "ddloops",3,... + % "trainloops",4,... + % "eq_parallelization_blocklength",eq_parallelization_blocklength, ... + % "eq_updatelatency",eq_updatelatency,... + % "eq_avg_blocklength",eq_avg_blocklength); + + [Eq_out] = eq.process(y_rx,y_digimod); + + results.EQ_out = Eq_out; + + % 3) digital demodulation object + digimod = PAMmapper(2^recorded_data.saveStructTemp.common.M,0); + + %estimated/ equalized symbol sequence + d_estimated = digimod.demap(Eq_out); + + %correct data symbols + d_correct = recorded_data.saveStructTemp.prms_out; + + % 4) BER calculation + % BER + [totalbits,errors,results.ber,loc] = calc_ber(d_estimated.signal(1:end,:),d_correct(1:end,:)',"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % [totalbits,errors,results.ber,loc] = calc_ber(d_estimated.signal(1:end,:) ,d_correct(1:end,:)',"skip",0,"returnErrorLocation",1); +end + +function std_per_lvl = calcleveldependentstd(rx_symbols,correct_symbols) + error_of_rx_signal = rx_symbols - correct_symbols; + levels = unique(correct_symbols); + for l = 1:4 + level_amplitude = levels(l); + std_per_lvl(l) = var(( 1/64 .* movsum(error_of_rx_signal(correct_symbols==level_amplitude),[64/2,64/2]) )); + %std_per_lvl(l) = std(error_of_rx_signal(correct_symbols==level_amplitude)); + end +end + +function scatterleveldependent(rx_symbols,correct_symbols,f_sym) + col = cbrewer2('Paired',8); + ccnt = -1; + figure1 = figure(); + levels = unique(correct_symbols); + + start = 1; + ende = length(correct_symbols); + + start = 30000; + ende = 40000; + + for l = 1:4 + ccnt = ccnt+2; + + level_amplitude = levels(l); + + symbols_for_lvl = NaN(1,length(correct_symbols)); + symbols_for_lvl(correct_symbols==level_amplitude) = rx_symbols(correct_symbols==level_amplitude); + std_lvl(l) = std(symbols_for_lvl,'omitnan'); + xax_in_sec = ((1:length(correct_symbols)) / f_sym) * 1e6; + xax_in_sec = 1:length(correct_symbols); + + scatter(xax_in_sec(start:ende),symbols_for_lvl(start:ende),10,'.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5,'MarkerEdgeColor',col(ccnt,:)); + hold on; + end + + std_lvl = round(std_lvl,2); + disp(std_lvl); + + ccnt = 0; + + % Add the windowed/ smoothed curves + for l = 1:4 + ccnt = ccnt+2; + level_amplitude = levels(l); + + symbols_for_lvl = NaN(1,length(correct_symbols)); + + movmean = 1/250 .* movsum(rx_symbols(correct_symbols==level_amplitude),[250/2,250/2]); + + symbols_for_lvl(correct_symbols==level_amplitude) = movmean; + + nanx = isnan(symbols_for_lvl); + t = 1:numel(symbols_for_lvl); + symbols_for_lvl(nanx) = interp1(t(~nanx), symbols_for_lvl(~nanx), t(nanx)); + + xax_in_sec = ((1:length(correct_symbols)) / f_sym) * 1e6; + xax_in_sec = 1:length(correct_symbols); + + plot(xax_in_sec(start:ende),symbols_for_lvl(start:ende),'Color',col(ccnt,:)); + hold on + end + + %yline(max(rx_symbols(correct_symbols==levels(2)))) + + if 0 + annotation(figure1,'textbox',... + [0.660523809523809 0.844444444444448 0.133523809523809 0.0603174603174607],... + 'String',['\sigma = ',num2str(std_lvl(4))],... + 'LineWidth',1.8,... + 'LineStyle','none',... + 'FontSize',12,... + 'FitBoxToText','off'); + + % Create textbox + annotation(figure1,'textbox',... + [0.667666666666665 0.642857142857147 0.133523809523809 0.0603174603174607],... + 'String',['\sigma = ',num2str(std_lvl(3))],... + 'LineWidth',1.8,... + 'LineStyle','none',... + 'FontSize',12,... + 'FitBoxToText','off'); + + % Create textbox + annotation(figure1,'textbox',... + [0.671238095238093 0.442857142857148 0.133523809523809 0.0603174603174608],... + 'String',['\sigma = ',num2str(std_lvl(2))],... + 'LineWidth',1.8,... + 'LineStyle','none',... + 'FontSize',12,... + 'FitBoxToText','off'); + + % Create textbox + annotation(figure1,'textbox',... + [0.670047619047616 0.265079365079371 0.133523809523809 0.0603174603174608],... + 'String',['\sigma = ',num2str(std_lvl(1))],... + 'LineWidth',1.8,... + 'LineStyle','none',... + 'FontSize',12,... + 'FitBoxToText','off'); + end + + xlim([0, 2.6]) + ylim([-2 2]) + xlabel('Time in $\mu$s'); + ylabel('Normalized Amplitude'); + +end + + diff --git a/projects/MPI_analysis_Mai24/simulation_entry.m b/projects/MPI_analysis_Mai24/simulation_entry.m new file mode 100644 index 0000000..7026634 --- /dev/null +++ b/projects/MPI_analysis_Mai24/simulation_entry.m @@ -0,0 +1,14 @@ + + +% settings + +M = 4; +datarate = 224e9; +kover = 8; + + +% construct trasnmit signal + +% transmit including MPI + +% receive \ No newline at end of file diff --git a/projects/MPI_offline_analysis/mpi_analysis_pam4.m b/projects/MPI_offline_analysis/mpi_analysis_pam4.m index a75b9f6..27e3158 100644 --- a/projects/MPI_offline_analysis/mpi_analysis_pam4.m +++ b/projects/MPI_offline_analysis/mpi_analysis_pam4.m @@ -11,10 +11,9 @@ end % SETTINGS - optimize_mudc = 0; -run_sir_sweep = 0; -run_feed_forward = 1; +run_sir_sweep = 1; +run_feed_forward = 0; run_baseline = 0; run_ideal_dc_tap = 0; plot_timesignal = 1; @@ -33,24 +32,26 @@ for bl = 1:numel(block_loop) mudc_loop = [0, 0.0001,0.0005, 0.001,0.005, 0.01:0.01:0.1, 0.2:0.1:1]; current_filename = 'pam4__loop_14_92Gbd_19092023_1513.mat'; - current_filename = 'pam4__loop_29_56Gbd_19092023_1546.mat'; + % current_filename = 'pam4__loop_29_56Gbd_19092023_1546.mat'; recorded_data = load([foldername,filesep, current_filename]); [best_mudc,ber]= optimizeMuDc(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc_loop); -% figure(16) -% hold on -% plot(mudc_loop,ber); -% set(gca, 'YScale', 'log'); -% set(gca, 'XScale', 'log'); -% scatter(best_mudc,ber(mudc_loop==best_mudc),100,'Marker','x','LineWidth',2) -% yline(ber(1)); -% xlim([0,1]) + figure(16) + hold on + scatter(mudc_loop(mudc_loop~=best_mudc),ber(mudc_loop~=best_mudc),'Marker','*','LineWidth',1); + set(gca, 'YScale', 'log'); + set(gca, 'XScale', 'log'); + scatter(best_mudc,ber(mudc_loop==best_mudc),100,'Marker','*','LineWidth',2) + yline(ber(1)); + xlim([0,1]) end if run_sir_sweep - - mudc = best_mudc; + + %mudc = 0.04;%best_mudc; + mudc = 0; + [ber,sir,fsym] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc); plotSirSweep(ber,sir,fsym,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc) @@ -67,7 +68,7 @@ end if run_feed_forward - current_filename = 'pam4__loop_14_92Gbd_v19092023_1513.mat'; + current_filename = 'pam4__loop_14_92Gbd_19092023_1513.mat'; recorded_data = load([foldername,filesep, current_filename]); eq_avg_blocklength = [50,100,1000,3000]; [best_block,ber]= optimizeAvgBlocklength(recorded_data,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,0); @@ -91,7 +92,8 @@ end if run_baseline - + + % baseline means no DC removal alg! mudc = 0; [ber_base,sir_base,fsym_base] = runSIRsweep(eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength,mudc); @@ -219,7 +221,7 @@ function [ber,sir,fsym] = runSIRsweep(eq_parallelization_blocklength,eq_updatela if ismac foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km'; else - foldername = 'C:\Users\Silas\Documents\MATLAB\Labor_Datensatz_PAM4_MPI\pam4_10km_1km'; + foldername = 'C:\Users\Silas\Documents\MATLAB\Datensätze\Labor_Datensatz_PAM4_MPI_OFC2023\pam4_10km_1km'; end allfiles = dir(foldername); @@ -236,7 +238,7 @@ function [ber,sir,fsym] = runSIRsweep(eq_parallelization_blocklength,eq_updatela results = runPostProcessing(recorded_data,mudc,eq_parallelization_blocklength,eq_updatelatency,eq_avg_blocklength); ber(i) = results.ber; fsym(i) = results.fsym; - sir(i) = results.sir + sir(i) = results.sir; disp(['fsym: ',num2str(results.fsym*1e-9),'GBd, SIR: ',num2str(results.sir),' ->> BER: ',sprintf('%2E',results.ber)]); @@ -293,13 +295,27 @@ function results = runPostProcessing(recorded_data,mudc,eq_parallelization_block y_rx.fs = 2.*results.fsym; % 0) build tx reference Signal for eq training - y_digimod = Electricalsignal(recorded_data.saveStructTemp.digi_mod_out'); + y_digimod = Informationsignal(recorded_data.saveStructTemp.digi_mod_out'); y_digimod.fs = results.fsym; % 1) normlaize y_rx = y_rx.normalize("mode","rms"); + % eq = EQ_silas_ofc("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... + % "sps",2,... + % "mu_dc_dd",mudc,... + % "mu_dc_train",mudc,... + % "mu_ffe_train",0,... + % "mu_dfe_train",0.005,... + % "mu_ffe_dd",[0.0004 0.0004 0.0004],... + % "mu_dfe_dd",0.005,... + % "ddloops",3,... + % "trainloops",4,... + % "eq_parallelization_blocklength",eq_parallelization_blocklength, ... + % "eq_updatelatency",eq_updatelatency,... + % "eq_avg_blocklength",eq_avg_blocklength); + eq = EQ_silas("Ne",[25,0,0],"Nb",[2,0,0],"trainlength",4096,... "sps",2,... "mu_dc_dd",mudc,... @@ -328,7 +344,9 @@ function results = runPostProcessing(recorded_data,mudc,eq_parallelization_block d_correct = recorded_data.saveStructTemp.prms_out; % 4) BER calculation - [totalbits,errors,results.ber,loc] = calc_ber(d_estimated.signal(1:end,:) ,d_correct(1:end,:)',"skip",0,"returnErrorLocation",1); + % BER + [totalbits,errors,results.ber,loc] = calc_ber(d_estimated.signal(1:end,:),d_correct(1:end,:)',"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % [totalbits,errors,results.ber,loc] = calc_ber(d_estimated.signal(1:end,:) ,d_correct(1:end,:)',"skip",0,"returnErrorLocation",1); end function std_per_lvl = calcleveldependentstd(rx_symbols,correct_symbols) diff --git a/projects/standard_system/imdd_minimal.m b/projects/standard_system/imdd_minimal.m new file mode 100644 index 0000000..1d51090 --- /dev/null +++ b/projects/standard_system/imdd_minimal.m @@ -0,0 +1,86 @@ + + +clear +M = 4; +fsym = 112e9; +fdac = 256e9; +kover = 8; +oneway_interference_meter = 0; +link_total_meter = 10000; +sir = 100; +rop = 0; + + +[D,B] = PAMsource("order",18,"useprbs",1,"fsym",fsym,"M",M).process(); + +if 1 + S = Pulseformer("fsym",fsym,"fdac",256e9,"pulse","rrc","pulselength",16,"rrcalpha",0.1).process(D); +end + +if 1 + min_ = min(D.signal) * 1.3 ; + max_ = max(D.signal) * 1.3 ; + S.signal = clip(S.signal,min_,max_); +end + +X = M8199A("kover",kover).process(S).*0.7222; + +u_pi = 2.9; +vbias = -0.85*u_pi; +[O,extmodlaser] = EML("mode",eml_mode.im_cosinus,"power",0,"fsimu",X.fs,"lambda",1310,"bias",vbias,"u_pi",u_pi,"linewidth",0,"randomkey",1).process(X); + +O.eye(fsym,M); + +O.spectrum("fignum",112,"displayname",'Opt. transmit spectrum'); + +if oneway_interference_meter ~= 0 + %%%%% Ping Pong/ Interference Path %%%%%% + I = Fiber("fsimu",O.fs,"fiber_length",oneway_interference_meter*2/1000,"alpha",0.3,"D",0,"lambda0",1320,"gamma",0,"Dslope",0.07).process(O); + I = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",-sir).process(I); + + %%%%% Delay the "main" signal as in reality the interference is "older" than the main signal %%%%%% + [O,dly] = O.delay("delay_meter",oneway_interference_meter*2); + + %%%%% Recombine Interference and Main Signal %%%%%% + O = O + I; + + %%%%% Cut (due to the delays there is a jump in the signals) %%%%%% + if dly == 0;dly = 1;end + O.signal = O.signal(ceil(dly):end); +end + +%%%%% Propagate through fiber %%%%%% +O = Fiber("fsimu",O.fs,"fiber_length",link_total_meter/1000,"alpha",0.3,"D",0,"lambda0",1320,"gamma",0,"Dslope",0.07).process(O); + +% Set ROP +O = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(O); + +E = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20).process(O); + +E = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(E); + +Scpe_sig = Scope("fsimu",fdac*kover,"fadc",256e9,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Filter('filtdegree',4,"f_cutoff",63e9,"fs",256e9,"filterType",filtertypes.butterworth,"active",true)).process(E); + +%%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",Scpe_sig.fs,"fs_out",2*fsym); + % Scpe_sig.plot('fignum',12345,'displayname','bla') + +%%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,delayed,cuts] = Scpe_sig.tsynch("reference",D,"fs_ref",fsym); + +[EQ_sig] = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",80,"sps",2,"decide",1).process(Scpe_sig,D); + +[~,errors_bm,ber,errors] = calc_ber(EQ_sig.signal,B.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + +disp(['BER: ',sprintf('%.1E',ber),' - - ROP: ',num2str(O.power),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); + + + + + + + + diff --git a/wh.mat b/wh.mat new file mode 100644 index 0000000..320d18a Binary files /dev/null and b/wh.mat differ