diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index a2f70ae..c2a7244 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -344,14 +344,14 @@ classdef Signal % 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 [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); w=w.*1e-9; else [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,"centered","power","mean"); - p_lin = smooth(p_lin,0.05,'rloess'); + % p_lin = smooth(p_lin,0.05,'rloess'); end @@ -376,10 +376,10 @@ classdef Signal if options.normalizeToNyquist==0 xlabel("Frequency in GHz"); - %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)); 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 xlabel("Normalized Frequency"); xlim([-pi, pi]); diff --git a/Classes/01_transmit/ChannelFreqResp.m b/Classes/01_transmit/ChannelFreqResp.m index e77cf43..f204112 100644 --- a/Classes/01_transmit/ChannelFreqResp.m +++ b/Classes/01_transmit/ChannelFreqResp.m @@ -241,10 +241,14 @@ classdef ChannelFreqResp < handle xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on; %%% plot for publication - figure(98989);hold all;box on;title('Magnitude Freq. Response'); - xlim([0.2 .5*max(obj.faxis)*1e-9]); - ylim([-20, 2]); - plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2); + figure(1234);hold all;box on;title('Magnitude Freq. Response'); + %xlim([0.2 .5*max(obj.faxis)*1e-9]); + ylim([-40, 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; diff --git a/projects/400G_FTN_setups/imdd_mpi_dsp.m b/projects/400G_FTN_setups/imdd_mpi_dsp.m index d9464c2..2016158 100644 --- a/projects/400G_FTN_setups/imdd_mpi_dsp.m +++ b/projects/400G_FTN_setups/imdd_mpi_dsp.m @@ -5,17 +5,17 @@ params = struct; params.M = [4]; params.datarate = [448]; params.rop = [0]; -params.sir = 40;%15:1:40; +params.sir = 45; params.random_key_laser_phase = 10:20; 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_channelapproach = 1; +db_channelapproach = 0; -laser_linewidth = 50e5; +laser_linewidth = 5e5; random_key_sequence = 15; random_key_laser_phase = 66; sir = 20; @@ -61,14 +61,14 @@ for M = wh.parameter.M.values % MAIN SIGNAL %%%%% 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,... "applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence,... "db_precode",db_precode,"db_encode",db_encode,... "mrds_code",usemrds,"mrds_blocklength",512).process(); %%%%% 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,... "applypulseform",0,"pulseformer",Pform,"randkey",random_key_sequence+1,... "db_precode",db_precode,"db_encode",db_encode,... @@ -185,8 +185,8 @@ for M = wh.parameter.M.values end % 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 %%%%%% 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 % ref symbols and transm. sequence are precoded [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 = Duobinary().decode(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); - + + + 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 [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); @@ -227,21 +233,39 @@ for M = wh.parameter.M.values elseif postfilter [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); - - % EQ_sig.plot("displayname",'After VNLE','fignum',90,'clear',1); - % Quantization is too far from orig. symbols -> - % error psd is quite different - % Sym_ = PAMmapper(M,0).quantize(EQ_sig); - % Noi_ = Sym_-EQ_sig; - % Noi_.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:)); + %%% REMOVE DC peak from Noi PSD + S = Noi.signal; + N1 = 1001; + + % 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_(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); [~,~,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,:)); - Noi.normalize('mode','rms').spectrum('displayname','Noise PSD optimal','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:)); + EQ_sig.spectrum('displayname','EQ Out','fignum',12345,'normalizeTo0dB',0,'normalizeToNyquist',0,'color',cols(1,:)); + + + Noi.spectrum('displayname','Noise PSD optimal','fignum',22,'normalizeTo0dB',1,'normalizeToNyquist',0,'color',cols(2,:)); for nc = 1:3 @@ -258,24 +282,16 @@ for M = wh.parameter.M.values if 1 cols = linspecer(12); + % EQ_sig_filt.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',0,'color',cols(nc+2,:)); - - EQ_sig_filt.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+2,:)); - - % [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs); - % h = h/max(abs(h)); - % hold on - % w_ = (w - Noi.fs/2); - % figure(123) - % plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); - - [h,w] = freqz(1,burg_coeff,length(Noi),"whole"); + [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs); + % h = 1./h; h = h/max(abs(h)); hold on - w_ = (w - pi); - plot(w_,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); + w_ = (w - Noi.fs/2); + figure(22) + plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']); - end Rx_bits = PAMmapper(M,0).demap(EQ_sig_mlse); @@ -286,41 +302,7 @@ for M = wh.parameter.M.values 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); @@ -379,16 +361,16 @@ ber_mlse=squeeze(mean(ber_mlse,1)); ber_vnle = mean(ber_vnle,1); % Create the initial plot -figure(44); +figure(466); a = gca; hold on; % Retain the plot so new points can be added without complete redraw dispname = ['Lw: ',num2str(laser_linewidth.*1e-6),' MHz']; - -plot(wh.parameter.sir.values,ber_vnle,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['VNLE ',dispname]); -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(:,2),"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['MLSE 2',dispname]); -plot(wh.parameter.sir.values,ber_mlse(:,3),"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['MLSE 3',dispname]); +cols = linspecer(6); +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",['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",['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",['PAM',num2str(M),'MLSE 3',dispname],'Color',cols(4,:)); yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); xlabel('Received Optical Power (dBm)'); diff --git a/projects/Lab_2024/offline_dsp_analysis/load_n_dsp.m b/projects/Lab_2024/offline_dsp_analysis/load_n_dsp.m index fe8434f..784127b 100644 --- a/projects/Lab_2024/offline_dsp_analysis/load_n_dsp.m +++ b/projects/Lab_2024/offline_dsp_analysis/load_n_dsp.m @@ -1,13 +1,14 @@ % load data points 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_v2'; 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'; end + stuff = load([foldername,filesep,filename,'wh']); %stuff = load([foldername,filesep,'PAM4_v2_10km_wh']); wh = stuff.obj; diff --git a/projects/Lab_2024/offline_dsp_analysis/load_n_plot_transfer_characteristic.m b/projects/Lab_2024/offline_dsp_analysis/load_n_plot_transfer_characteristic.m new file mode 100644 index 0000000..a23d25a --- /dev/null +++ b/projects/Lab_2024/offline_dsp_analysis/load_n_plot_transfer_characteristic.m @@ -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();