diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 0f3578e..5905556 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -757,11 +757,11 @@ classdef Signal pkpos = sort(pkpos); - if mean(w) > 10 || mean(p) > 15 - return - else - sequenceFound = 1; - end + % if mean(w) > 15 || mean(p) > 15 + % return + % else + % sequenceFound = 1; + % end if options.debug_plots figure(121212);clf diff --git a/Classes/01_transmit/Pulseformer.m b/Classes/01_transmit/Pulseformer.m index 4efe336..d137924 100644 --- a/Classes/01_transmit/Pulseformer.m +++ b/Classes/01_transmit/Pulseformer.m @@ -43,7 +43,7 @@ classdef Pulseformer 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.signal = obj.process_(signalclass_in); % append to logbook lbdesc = 'Applied Pulseshaping'; @@ -65,109 +65,171 @@ classdef Pulseformer % Cant be seen from outside! So put all your functions here that can/ % shall not be called from outside - function data_out = process_(obj,data_in) - %METHOD1 Summary of this method goes here - % Detailed explanation goes here - arguments(Input) - obj - data_in - end - - arguments(Output) - data_out - end - - if ~rem(obj.fdac,obj.fsym) - %ist ein Vielfaches - sps = obj.fdac / obj.fsym; - p = sps; - q = 1; - else - %ist kein Vielfaches - p = obj.fsym / gcd(obj.fdac, obj.fsym); %upsampling p->->-> - q = obj.fdac/ gcd(obj.fdac, obj.fsym); %downsampling <-q - sps= q; %sps während dem pulse shaping - end - - if obj.pulse == pulseform.rc - filtertype = 'normal'; - elseif pulseform.rrc - filtertype = 'sqrt'; - end - - %Bau das Filter (hier rc) - racos_len = obj.pulselength*2; - h = rcosdesign(obj.alpha,racos_len,sps,filtertype); - % h = h./ max(h); - - if obj.matched - h = conj(fliplr(h)); - end - - manual_cyclic_convolution = 0; - upfirdn_convolution = 1; - - if manual_cyclic_convolution - - % Apply filter the long way (from move_it) - data_in = data_in'; - blen = length(data_in)*sps; - - % oversample symbol sequence - symbolov=zeros(size(data_in,1),blen); - symbolov(:,1:sps:blen-sps+1)=data_in; - H=fft(h,blen); - - % Convolution of Bit sequence with impulse response - data_out=ifft( fft(symbolov.') .* repmat( H,size(data_in,1),1 ).' ).'; - data_out = circshift(data_out,[0 -(obj.pulselength*sps)]); +function data_out = process_(obj, data_in_signal) + % Extract the incoming sampling rate + % Safety check: If fs is missing (e.g. raw symbols), assume it is fsym + if isprop(data_in_signal, 'fs') && ~isempty(data_in_signal.fs) + f_in = data_in_signal.fs; + else + f_in = obj.fsym; + end - if rem(obj.fdac,obj.fsym) - data_out = data_out(1:q:end); - end + f_out = obj.fdac; - end + % 1. Calculate Resampling Factors (P and Q) + % We need rational approximation: f_out/f_in = p/q + [p, q] = rat(f_out / f_in); - if upfirdn_convolution + % 2. Calculate SPS for the Filter Design + % The filter operates at the INTERMEDIATE rate (f_in * p). + % We need to know how many samples represent one symbol AT THAT RATE. + fs_intermediate = f_in * p; + sps_filter = fs_intermediate / obj.fsym; - %Apply Filter using Matlab build in fctn. + % 3. Filter Design + if obj.pulse == pulseform.rc + filtertype = 'normal'; + elseif obj.pulse == pulseform.rrc % assuming enum logic holds + filtertype = 'sqrt'; + end - data_out_ = upfirdn(data_in,h,p,q); + % Standard RRC Design + % Note: rcosdesign sps must be integer? Usually yes, but for polyphase it can handle it. + % If sps_filter is not integer, rcosdesign might complain. + % For your setup (powers of 2), it will likely be integer. + racos_len = obj.pulselength; % span in symbols + h = rcosdesign(obj.alpha, racos_len, sps_filter, filtertype); - %cut signal, which is longer due to fir filter - st = round(p/q*racos_len/2); %we need to cut y_out - en = round(st + (length(data_in)*p/q)); - data_out = data_out_(st:en); + % Matched Filter Flip (Complex Conjugate Time Reversal) + if obj.matched + h = conj(fliplr(h)); + end - end + % 4. Processing + % Apply upfirdn using the calculated P and Q + data_out_ = upfirdn(data_in_signal.signal, h, p, q); - if upfirdn_convolution && manual_cyclic_convolution - figure() - subplot(2,1,1) - title("Convolution vs. Upfirdn and Cut") - hold on - % plot(data_out_(1:200),'DisplayName','Matlab upfirdn'); - plot(data_out(1:200),'DisplayName','By Hand cyclic convolution') - subplot(2,1,2) - hold on - plot(data_out(1:2000),'DisplayName','OUT'); - plot(data_in(1:2000),'DisplayName','IN'); - end + % 5. Trim Tail (Group Delay Correction) + % The delay of linear phase filter is (N-1)/2 samples @ intermediate rate + delay_samples_intermediate = (length(h) - 1) / 2; + + % Convert delay to output samples + delay_samples_out = delay_samples_intermediate / q; + + % We usually want to trim the "start" transient + st = floor(delay_samples_out) + 1; + + % Calculate expected output length + len_out = ceil(length(data_in_signal.signal) * p / q); + + % Cut + data_out = data_out_(st : st + len_out - 1); +end - %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 -% data_out = data_out./scale; - - % data_out = data_out'; - - %Check output integrity - if abs(round(p/q * length(data_in)) - length(data_out)) > 4 - warning('Check signal length after pulse shaping'); - %disp('Check signal length after pulse shaping'); - end - - - end +% +% function data_out = process_(obj,data_in) +% %METHOD1 Summary of this method goes here +% % Detailed explanation goes here +% arguments(Input) +% obj +% data_in +% end +% +% arguments(Output) +% data_out +% end +% +% if ~rem(obj.fdac,obj.fsym) +% %ist ein Vielfaches +% sps = obj.fdac / obj.fsym; +% p = sps; +% q = 1; +% else +% %ist kein Vielfaches +% p = obj.fsym / gcd(obj.fdac, obj.fsym); %upsampling p->->-> +% q = obj.fdac/ gcd(obj.fdac, obj.fsym); %downsampling <-q +% sps= q; %sps während dem pulse shaping +% end +% +% if obj.pulse == pulseform.rc +% filtertype = 'normal'; +% elseif pulseform.rrc +% filtertype = 'sqrt'; +% end +% +% %Bau das Filter (hier rc) +% racos_len = obj.pulselength*2; +% h = rcosdesign(obj.alpha,racos_len,sps,filtertype); +% % h = h./ max(h); +% +% if obj.matched +% h = conj(fliplr(h)); +% end +% +% manual_cyclic_convolution = 0; +% upfirdn_convolution = 1; +% +% if manual_cyclic_convolution +% +% % Apply filter the long way (from move_it) +% data_in = data_in'; +% blen = length(data_in)*sps; +% +% % oversample symbol sequence +% symbolov=zeros(size(data_in,1),blen); +% symbolov(:,1:sps:blen-sps+1)=data_in; +% H=fft(h,blen); +% +% % Convolution of Bit sequence with impulse response +% data_out=ifft( fft(symbolov.') .* repmat( H,size(data_in,1),1 ).' ).'; +% data_out = circshift(data_out,[0 -(obj.pulselength*sps)]); +% +% if rem(obj.fdac,obj.fsym) +% data_out = data_out(1:q:end); +% end +% +% end +% +% if upfirdn_convolution +% +% %Apply Filter using Matlab build in fctn. +% +% data_out_ = upfirdn(data_in,h,p,q); +% +% %cut signal, which is longer due to fir filter +% st = round(p/q*racos_len/2); %we need to cut y_out +% en = round(st + (length(data_in)*p/q)); +% data_out = data_out_(st:en); +% +% end +% +% if upfirdn_convolution && manual_cyclic_convolution +% figure() +% subplot(2,1,1) +% title("Convolution vs. Upfirdn and Cut") +% hold on +% % plot(data_out_(1:200),'DisplayName','Matlab upfirdn'); +% plot(data_out(1:200),'DisplayName','By Hand cyclic convolution') +% subplot(2,1,2) +% hold on +% plot(data_out(1:2000),'DisplayName','OUT'); +% plot(data_in(1:2000),'DisplayName','IN'); +% end +% +% %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 +% % data_out = data_out./scale; +% +% % data_out = data_out'; +% +% %Check output integrity +% if abs(round(p/q * length(data_in)) - length(data_out)) > 4 +% warning('Check signal length after pulse shaping'); +% %disp('Check signal length after pulse shaping'); +% end +% +% +% end end diff --git a/Classes/Moveit_wrapper.m b/Classes/Moveit_wrapper.m index 621ab67..8c7ea8a 100644 --- a/Classes/Moveit_wrapper.m +++ b/Classes/Moveit_wrapper.m @@ -27,7 +27,7 @@ classdef Moveit_wrapper < handle % Ensure the function exists if ~exist(obj.moveit_function_name, 'file') - error('Function "%s" does not exist.', obj.moveit_function_name); + error('Function "%s" does not exist. The move-it module must be on path for Matlab, otherwise the wrapper can not call it...', obj.moveit_function_name); end % Step 1: Get default parameters by calling moveit module diff --git a/Functions/EQ_structures/ffe.m b/Functions/EQ_structures/ffe.m index 6cfefc5..df4138e 100644 --- a/Functions/EQ_structures/ffe.m +++ b/Functions/EQ_structures/ffe.m @@ -131,11 +131,11 @@ switch precode_mode tx_bits_precoded = mapper.demap(tx_symbols_precoded); rx_bits = mapper.demap(eq_signal_hd_precoded); - [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); % B) Just determine BER rx_bits = mapper.demap(eq_signal_hd); - [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); case db_mode.db_precoded % Data is precoded on TX side @@ -143,12 +143,12 @@ switch precode_mode eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd, "M", M); eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M); rx_bits_decoded = mapper.demap(eq_signal_hd_decoded); - [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); % B) Omit the Coding by comparing with demapped TX symbol sequence tx_bits_demapped = mapper.demap(tx_symbols); rx_bits = mapper.demap(eq_signal_hd); - [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); end end diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m index 530e692..c535c19 100644 --- a/projects/FSO_transmission/first_analysis.m +++ b/projects/FSO_transmission/first_analysis.m @@ -1,12 +1,141 @@ -datas = load("C:\Users\Silas\Downloads\FSO_FP_QCL_60umUTC\FSO_FP_QCL_60umUTC\6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=225mA_RoP=28.02mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +base = "C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC"; +mode = 0; %0 oder 1 +M = 2; -fsym = 6e9; -fs = 80e9; +all_files = dir(fullfile(base, "**/*.mat")); -rawsignal = datas.tr.lastData(1).trace.ch3.RawData; +if M == 2 + tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); + filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); +elseif M == 4 + tx_data = load("C:\Users\Silas\Nextcloud\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); + filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); +end -Scope_sig = Electricalsignal(rawsignal,"fs",fs); +if mode == 1 + [f, p] = uigetfile(fullfile(base, "**/*.mat")); + if f~=0 + filename = fullfile(p,f); + end +end + + +datas = load(filename); + +%% +str = filename; +M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); +assert(M==M_); +fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); +fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); +I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); +rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); +L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); +pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); +rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); +mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); + +%% +% Tx data + +Bits = Informationsignal(tx_data.tx_data,"fs",fsym); +Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); + +mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there +PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme + +Symbols_ = PM.map(Bits) .* PM.scaling; +assert(isequal(Symbols.signal,Symbols_.signal)); + +Bits_ = PM.demap(Symbols); +[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); +assert(ber == 0); + +%% For comparison, apply pulsef on Tx Symbols +Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rc","pulselength",16,"alpha",rolloff); +Digi_sig_compare = Pform.process(Symbols); +MF = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); +Rx_sig_compare = MF.process(Digi_sig_compare); + +%% + +% Rx Data +traceData = datas.tr.lastData(2).trace.ch3; + +%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin +scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; +demystified = isequal(traceData.YData,scoperead_volts); +assert(demystified); + +Scope_sig = Electricalsignal(traceData.YData,"fs",fs); Scope_sig.plot("displayname",'raw','fignum',100); +Scope_sig.spectrum("displayname",'raw','fignum',101) +% 1) matched filter +% pulse is symmetric, hence we can use pulsef firectly as matched filter. +% It feels off (bit I think correct) that the fsym is now the output freq.!! +% -> output 2 sps to omit timing recovery!? +Pform = Pulseformer("fsym",fsym,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); +Rx_matched = Pform.process(Scope_sig); +Rx_matched.spectrum("displayname",'Signal after matched filter','fignum',1); + +% timing sync -> at this point we still have no symbol timing recovery, we +% try to do this with 2sps EQ! +[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + +%% not working.. +Rx_synced = Rx_synced_cell{1}; +len_tr = 4096*2; +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; +duob_mode = db_mode.no_db; + +Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',2); + +Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',0); +Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',0); +Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',0); + +if M == 2 + ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature +elseif M == 4 + ber_in_paper = 10^(-2.5); +end + +%% -------------------- FFE -------------------- +% requires some more digging what is going on :-) +eq_ffe = EQ("Ne",[50, 5, 5],"Nb",[2,0,0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); + +ffe_results = ffe(eq_ffe,M,Rx_synced,Symbols,Bits, ... + "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ... + "eth_style_symbol_mapping",mapping_style); + +% ffe_results.metrics.print +fprintf('My EQ: %.1e \n',ffe_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); + +%% -------------------- VNLE + MLSE -------------------- + +pf_ncoeffs = 1; +eq_v = EQ("Ne",[100, 5, 5],"Nb",[0, 0, 0], ... + "training_length",len_tr,"training_loops",5,"dd_loops",5, ... + "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... + "FFEmu",0,"plotfinal",0,"ideal_dfe",1); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); + +[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, M, Rx_synced, Symbols, Bits, ... + "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); + +mlse_results.metrics.print +fprintf('My EQ: %.1e \n',mlse_results.metrics.BER); +fprintf('Paper: %.1e \n \n',ber_in_paper); \ No newline at end of file diff --git a/test/matched_filter_minimal.m b/test/matched_filter_minimal.m index b298fac..3f2ff9f 100644 --- a/test/matched_filter_minimal.m +++ b/test/matched_filter_minimal.m @@ -13,7 +13,6 @@ db_precode = 0; db_encode = 0; -rcalpha = 0.05; kover = 16; vbias_rel = 0.5; u_pi = 2.9; @@ -47,15 +46,21 @@ if M == 6 bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); end -bits = Informationsignal(bitpattern); +bits = Informationsignal(bitpattern.'); symbols = PAMmapper(M,0).map(bits); symbols.fs = fsym; -Pform = Pulseformer("fsym",fsym,"fs",fsym,"alpha",0.6,"pulse","rrc","pulselength",64,"matched_sps",4,"output_sps",2); -Pform.process(symbols); +symbols.spectrum("displayname",'Symbols','fignum',1); +%% RRC Shaping +rcalpha = 1; +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); +Digi_sig = Pform.process(symbols); +Digi_sig.spectrum("displayname",'Signal after pluse shaping','fignum',1); + +%% RRC Matched Filtering +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); +Rx_sig = Pform.process(Digi_sig); +Rx_sig.spectrum("displayname",'Signal after matched filter','fignum',1); -MF = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",1024,"alpha",rcalpha,"matched",0); -Digi_sig = MF.process(Digi_sig); -Digi_sig.spectrum("displayname",'Signal after shaping','fignum',1);