CHnages from silas PC

This commit is contained in:
sioe
2024-11-04 09:00:51 +01:00
parent 6dca327d10
commit 9ed763d094
5 changed files with 80 additions and 85 deletions

View File

@@ -344,14 +344,14 @@ classdef Signal
% spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname); % spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname);
N = 2^(nextpow2(length(obj.signal))-10); N = 2^(nextpow2(length(obj.signal))-2);
if options.normalizeToNyquist==0 if options.normalizeToNyquist==0
[p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean");
w=w.*1e-9; w=w.*1e-9;
else else
[p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,"centered","power","mean"); [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,"centered","power","mean");
p_lin = smooth(p_lin,0.05,'rloess'); % p_lin = smooth(p_lin,0.05,'rloess');
end end
@@ -379,7 +379,7 @@ classdef Signal
% xlim([-obj.fs/2 obj.fs/2].*1e-9) % xlim([-obj.fs/2 obj.fs/2].*1e-9)
edgetick = 2^(nextpow2(obj.fs*1e-9)); edgetick = 2^(nextpow2(obj.fs*1e-9));
xticks([-edgetick:16:edgetick]); xticks([-edgetick:16:edgetick]);
xlim([100*round( min(w.*1e-9)/100,1)-10,100*round( max(w.*1e-9)/100,1)+10]) xlim([100*round( min(w)/100,1)-10,100*round( max(w)/100,1)+10])
else else
xlabel("Normalized Frequency"); xlabel("Normalized Frequency");
xlim([-pi, pi]); xlim([-pi, pi]);

View File

@@ -241,10 +241,14 @@ classdef ChannelFreqResp < handle
xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on; xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;
%%% plot for publication %%% plot for publication
figure(98989);hold all;box on;title('Magnitude Freq. Response'); figure(1234);hold all;box on;title('Magnitude Freq. Response');
xlim([0.2 .5*max(obj.faxis)*1e-9]); %xlim([0.2 .5*max(obj.faxis)*1e-9]);
ylim([-20, 2]); ylim([-40, 2]);
plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2); Havg_smooth = smooth(Havg,50);
symaxis = (obj.faxis-(obj.f_ref/2))/1e9;
Havg = fftshift(Havg);
Havg = smooth(Havg);
plot(symaxis, 20*log10(abs(Havg)),'LineWidth',0.5);
grid on; grid on;

View File

@@ -5,17 +5,17 @@ params = struct;
params.M = [4]; params.M = [4];
params.datarate = [448]; params.datarate = [448];
params.rop = [0]; params.rop = [0];
params.sir = 40;%15:1:40; params.sir = 45;
params.random_key_laser_phase = 10:20; params.random_key_laser_phase = 10:20;
precomp_mode = 0; %0=do nothing ; 1= measure; 2=precomp active precomp_mode = 0; %0=do nothing ; 1= measure; 2=precomp active
postfilter = 0; % noise whiten. approach -> Postfilter + MLSE postfilter = 1; % noise whiten. approach -> Postfilter + MLSE
db_precode = 1; db_precode = 0;
db_encode = 0; db_encode = 0;
db_channelapproach = 1; db_channelapproach = 0;
laser_linewidth = 50e5; laser_linewidth = 5e5;
random_key_sequence = 15; random_key_sequence = 15;
random_key_laser_phase = 66; random_key_laser_phase = 66;
sir = 20; sir = 20;
@@ -61,14 +61,14 @@ for M = wh.parameter.M.values
% MAIN SIGNAL % MAIN SIGNAL
%%%%% Symbol Generation MAIN %%%%%% %%%%% Symbol Generation MAIN %%%%%%
[Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",19,"useprbs",1,...
"fs_out",M8199.fdac,"applyclipping",0,"clipfactor",1.5,... "fs_out",M8199.fdac,"applyclipping",0,"clipfactor",1.5,...
"applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence,... "applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence,...
"db_precode",db_precode,"db_encode",db_encode,... "db_precode",db_precode,"db_encode",db_encode,...
"mrds_code",usemrds,"mrds_blocklength",512).process(); "mrds_code",usemrds,"mrds_blocklength",512).process();
%%%%% Symbol Generation INTERFERENCE %%%%%% %%%%% Symbol Generation INTERFERENCE %%%%%%
[Digi_sig_I,Symbols_I,Bits_I] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,... [Digi_sig_I,Symbols_I,Bits_I] = PAMsource("fsym",fsym,"M",M,"order",19,"useprbs",0,...
"fs_out",M8199.fdac,"applyclipping",0,"clipfactor",1.5,... "fs_out",M8199.fdac,"applyclipping",0,"clipfactor",1.5,...
"applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence+1,... "applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence+1,...
"db_precode",db_precode,"db_encode",db_encode,... "db_precode",db_precode,"db_encode",db_encode,...
@@ -185,8 +185,8 @@ for M = wh.parameter.M.values
end end
% Scpe_sig_normalized = Scpe_sig.normalize("mode","rms"); % Scpe_sig_normalized = Scpe_sig.normalize("mode","rms");
%
% Scpe_sig.normalize("mode","rms").spectrum("displayname",'After Scope','fignum',10); % Scpe_sig_normalized.normalize("mode","rms").spectrum("displayname",'After Scope','fignum',23);
%%%%%% Sample to 2x fsym %%%%%% %%%%%% Sample to 2x fsym %%%%%%
Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym);
@@ -213,11 +213,17 @@ for M = wh.parameter.M.values
if db_channelapproach if db_channelapproach
% ref symbols and transm. sequence are precoded % ref symbols and transm. sequence are precoded
[EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols));
EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig);
EQ_sig = Duobinary().decode(EQ_sig); EQ_sig = Duobinary().decode(EQ_sig);
Rx_bits = PAMmapper(M,0).demap(EQ_sig); Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,~,ber_vnle(i,j),~] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); [~,~,ber_vnle(i,j),~] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
EQ_sig.spectrum('displayname','EQ DB Out','fignum',12345,'normalizeTo0dB',0,'normalizeToNyquist',0,'color',cols(3,:));
Noi.spectrum('displayname','Noise PSD optimal','fignum',1234,'normalizeTo0dB',0,'normalizeToNyquist',0,'color',cols(4,:));
elseif db_encode elseif db_encode
[EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols);
@@ -228,20 +234,38 @@ for M = wh.parameter.M.values
[EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols);
% EQ_sig.plot("displayname",'After VNLE','fignum',90,'clear',1); %%% REMOVE DC peak from Noi PSD
S = Noi.signal;
N1 = 1001;
% Quantization is too far from orig. symbols -> % recursion
% error psd is quite different % Initialize the moving sum for the first window
% Sym_ = PAMmapper(M,0).quantize(EQ_sig); half_window = (N1 - 1) / 2;
% Noi_ = Sym_-EQ_sig; moving_sum = sum(S(1:N1));
% Noi_.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:));
% Calculate the first element of R1
S_(1:half_window+1) = S(1:half_window+1) - (moving_sum / N1);
% Loop over the signal and apply the recursive moving average subtraction
for n = (half_window+2):(length(S)-half_window)
% Update the moving sum by subtracting the oldest value and adding the new one
moving_sum = moving_sum - S(n-half_window-1) + S(n+half_window);
% Calculate the new value of R1
S_(n) = S(n) - (moving_sum / N1);
end
S_(n+1:length(S)) = S(n+1:end) - (moving_sum / N1);
Noi.signal = S_;
%%% END REMOVE DC PEAK %%%
Rx_bits = PAMmapper(M,0).demap(EQ_sig); Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,~,ber_vnle(i,j),~] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); [~,~,ber_vnle(i,j),~] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
EQ_sig.normalize('mode','rms').spectrum('displayname','EQ Out','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc,:)); EQ_sig.spectrum('displayname','EQ Out','fignum',12345,'normalizeTo0dB',0,'normalizeToNyquist',0,'color',cols(1,:));
Noi.normalize('mode','rms').spectrum('displayname','Noise PSD optimal','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:));
Noi.spectrum('displayname','Noise PSD optimal','fignum',22,'normalizeTo0dB',1,'normalizeToNyquist',0,'color',cols(2,:));
for nc = 1:3 for nc = 1:3
@@ -258,23 +282,15 @@ for M = wh.parameter.M.values
if 1 if 1
cols = linspecer(12); cols = linspecer(12);
% EQ_sig_filt.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',0,'color',cols(nc+2,:));
[h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs);
EQ_sig_filt.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+2,:)); % h = 1./h;
% [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs);
% h = h/max(abs(h));
% hold on
% w_ = (w - Noi.fs/2);
% figure(123)
% plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']);
[h,w] = freqz(1,burg_coeff,length(Noi),"whole");
h = h/max(abs(h)); h = h/max(abs(h));
hold on hold on
w_ = (w - pi); w_ = (w - Noi.fs/2);
plot(w_,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); figure(22)
plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']);
end end
@@ -286,41 +302,7 @@ for M = wh.parameter.M.values
else else
% S = Scpe_sig.signal;
% N1 = 101;
%
% % Initialize the running sum with the first window's sum
% running_sum = mean( S(1:N1) );
%
% % Calculate the first output value
% S_(1) = S(1) - running_sum;
%
% % Recursive running sum filter
% for n = 2 : length(S) - N1
% % Update running sum by removing the oldest sample and adding the newest
% avg_win(n) = mean( S(n:n+N1) );
% S_(n) = S(n) - avg_win(n);
% end
%
% % movmean
% S__ = S - movmean(S,[floor(N1/2),ceil(N1/2)]);
%
% % recursion
% % Initialize the moving sum for the first window
% half_window = (N1 - 1) / 2;
% moving_sum = sum(S(1:N1));
%
% % Calculate the first element of R1
% S___(half_window+1) = S(half_window+1) - (moving_sum / N1);
%
% % Loop over the signal and apply the recursive moving average subtraction
% for n = (half_window+2):(length(S)-half_window)
% % Update the moving sum by subtracting the oldest value and adding the new one
% moving_sum = moving_sum - S(n-half_window-1) + S(n+half_window);
%
% % Calculate the new value of R1
% S___(n) = S(n) - (moving_sum / N1);
% end
[EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols);
@@ -379,16 +361,16 @@ ber_mlse=squeeze(mean(ber_mlse,1));
ber_vnle = mean(ber_vnle,1); ber_vnle = mean(ber_vnle,1);
% Create the initial plot % Create the initial plot
figure(44); figure(466);
a = gca; a = gca;
hold on; % Retain the plot so new points can be added without complete redraw hold on; % Retain the plot so new points can be added without complete redraw
dispname = ['Lw: ',num2str(laser_linewidth.*1e-6),' MHz']; dispname = ['Lw: ',num2str(laser_linewidth.*1e-6),' MHz'];
cols = linspecer(6);
plot(wh.parameter.sir.values,ber_vnle,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['VNLE ',dispname]); plot(wh.parameter.sir.values,ber_vnle,"LineWidth",0.5,"LineStyle","--","Marker",".","MarkerSize",15,"DisplayName",['PAM',num2str(M),' VNLE ',dispname],'Color',cols(1,:));
plot(wh.parameter.sir.values,ber_mlse(:,1),"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['MLSE 1 ',dispname]); plot(wh.parameter.sir.values,ber_mlse(:,1),"LineWidth",0.5,"LineStyle","--","Marker",".","MarkerSize",15,"DisplayName",['PAM',num2str(M),'MLSE 1 ',dispname],'Color',cols(2,:));
plot(wh.parameter.sir.values,ber_mlse(:,2),"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['MLSE 2',dispname]); plot(wh.parameter.sir.values,ber_mlse(:,2),"LineWidth",0.5,"LineStyle","--","Marker",".","MarkerSize",15,"DisplayName",['PAM',num2str(M),'MLSE 2',dispname],'Color',cols(3,:));
plot(wh.parameter.sir.values,ber_mlse(:,3),"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['MLSE 3',dispname]); plot(wh.parameter.sir.values,ber_mlse(:,3),"LineWidth",0.5,"LineStyle","--","Marker",".","MarkerSize",15,"DisplayName",['PAM',num2str(M),'MLSE 3',dispname],'Color',cols(4,:));
yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off');
xlabel('Received Optical Power (dBm)'); xlabel('Received Optical Power (dBm)');

View File

@@ -1,13 +1,14 @@
% load data points % load data points
if 0 if 0
foldername = '/Users/silasoettinghaus/Nextcloud/Dokumente/02_Ablage_Office/Lab_Data_24/sir_sweep_sd40/'; foldername = 'C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\Lab_Data_24\sir_sweep_sd40';
filename = 'PAM4_56_'; filename = 'PAM4_56_';
filename = 'PAM4_56_v2'; filename = 'PAM4_56_v2';
else else
foldername = '/Users/silasoettinghaus/Nextcloud/Dokumente/02_Ablage_Office/Lab_Data_24/sir_sweep_sd40/DB_coded'; foldername = 'C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\Lab_Data_24\sir_sweep_sd40\DB_coded';
filename = 'PAM4_68_v2'; filename = 'PAM4_68_v2';
end end
stuff = load([foldername,filesep,filename,'wh']); stuff = load([foldername,filesep,filename,'wh']);
%stuff = load([foldername,filesep,'PAM4_v2_10km_wh']); %stuff = load([foldername,filesep,'PAM4_v2_10km_wh']);
wh = stuff.obj; wh = stuff.obj;

View File

@@ -0,0 +1,8 @@
precomp_path = "C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\standard_system";
precomp_filename = "lab_mpi_setup_2";
freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',92e9);
freqresp.load('loadPath',precomp_path,'fileName',precomp_filename);
freqresp.plot();