silas/imdd_silas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

changes from zurich; my mac

Browse Source
This commit is contained in:
Silas
2025-02-20 05:31:05 +01:00
parent d099efea03
commit cf94b67828
13 changed files with 197 additions and 148 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
Classes/00_signals/Signal.m
View File

@@ -355,14 +355,14 @@ classdef Signal
end
if isempty(options.fft_length)
options.fft_length = 2^(nextpow2(length(obj.signal))-7);
options.fft_length = 2^(nextpow2(length(obj.signal))-6);
end
if options.normalizeToNyquist == 0
[p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,obj.fs,"centered","psd","mean");
[p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,obj.fs,"centered","power","mean");
w = w.*1e-9;
else
[p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,"centered","psd","mean");
[p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,"centered","power","mean");
end
if options.normalizeTo0dB
@@ -411,8 +411,9 @@ classdef Signal
try
ylim([max(min(floor(min(p_dbm))-3, ax.YLim(1)),-40), min(max(ceil(max(p_dbm))+3, ax.YLim(2)),10)]);
catch
ylim([min(floor(min(p_dbm))-3, ax.YLim(1)), max(ceil(max(p_dbm))+3, ax.YLim(2))]);
ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]);
end
ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]);
yticks(-200:10:10);
grid on; grid minor;
legend
@@ -661,6 +662,7 @@ classdef Signal
obj Signal
options.reference Signal
options.fs_ref = 0;
options.debug_plots = 0;
end
%normalize the signal
@@ -680,6 +682,12 @@ classdef Signal
%estimate start pos of signal
maxpeaknum = floor(length(a)/length(b));
[pks,pkpos] = findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend');
if options.debug_plots
figure()
findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend')
end
shifts = lags(pkpos);
shifts = shifts(shifts>=0);
@@ -722,8 +730,7 @@ classdef Signal
end
%plot all synced signals and the ref signal
debug = 0;
if debug
if options.debug_plots
figure;hold on;
for i = 1:size(S,1)
plot(S{i}.normalize('mode','oneone').signal(1000:1100),'LineWidth',0.1,'Color',[0.2157 0.4941 0.7216]);
@@ -837,9 +844,9 @@ classdef Signal
mode = 1;
histpoints = 1024; %% verticale resolution
histpoints = 1024*2; %% 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
histpoints_horizontal = 512*2; %% horizontal resolution
hist_data=zeros(histpoints,histpoints_horizontal ); %% initilize eye diagram
if isa(obj,'Opticalsignal')
@@ -850,7 +857,10 @@ classdef Signal
sig = obj.signal;
end
x = (sig); %% make input signal rea)l
startpos = floor(0.1*length(sig));
endpos = floor(0.9*length(sig));
endpos = min(endpos,startpos+200000);
x = sig(startpos:endpos); %% make input signal rea)l
x = resample(x,fsym*histpoints_horizontal/2,obj.fs); %% up sample to original fsym rate

73
Classes/01_transmit/PAMmapper.m
View File

@@ -8,13 +8,21 @@ classdef PAMmapper
thresholds
levels
scaling
eth_style
end
methods
function obj = PAMmapper(M, unipolar)
function obj = PAMmapper(M, unipolar, options)
%PAMMAPPER Construct an instance of this class
% Detailed explanation goes here
arguments
M
unipolar
options.eth_style = 0;
end
obj.M = M;
obj.unipolar = unipolar;
obj.thresholds = obj.get_demodulation_thresholds();
@@ -23,6 +31,8 @@ classdef PAMmapper
obj.scaling = rms(obj.get_levels());
obj.eth_style = options.eth_style;
end
function out = map(obj,signal_in)
@@ -66,18 +76,25 @@ classdef PAMmapper
switch obj.M
case 2
% 2-ASK: BPSK / OOK
pam_sig=bitpattern(:,1);
if obj.unipolar==0
pam_sig=2*pam_sig-1;
if ~obj.eth_style
pam_sig = bitpattern(:,1);
if obj.unipolar==0
pam_sig=2*pam_sig-1;
end
else
pam_sig = -2*bitpattern(:,1) + 1;
end
case 4
% 4-ASK:
pam_sig=2*bitpattern(:,1)+(bitpattern(:,1)==bitpattern(:,2));
if obj.unipolar==0
pam_sig=2*pam_sig-3;
if ~obj.eth_style
pam_sig=2*bitpattern(:,1)+(bitpattern(:,1)==bitpattern(:,2));
if obj.unipolar==0
pam_sig=2*pam_sig-3;
end
else
pam_sig = (2*bitpattern(:,1)-1).*(-2*bitpattern(:,2)+3);
end
pam_sig = pam_sig/sqrt(5);
@@ -97,15 +114,22 @@ classdef PAMmapper
pam_sig = pam_sig/sqrt(10);
case 8
% 8-ASK:
x1 = bitpattern(:,1);
x2 = (bitpattern(:,1)==bitpattern(:,3));
x3 = x2~=bitpattern(:,2);
if ~obj.eth_style
x1 = bitpattern(:,1);
x2 = (bitpattern(:,1)==bitpattern(:,3));
x3 = x2~=bitpattern(:,2);
pam_sig = 4*x1 + 2*x2 + x3;
if obj.unipolar==0
pam_sig=2*pam_sig-7;
end
else
pam_sig = 4*x1 + 2*x2 + x3;
pam_sig = (bitpattern(:,1)*2-1).*(4+(2*bitpattern(:,2)-1).*(-2*bitpattern(:,3)+3));
if obj.unipolar==0
pam_sig=2*pam_sig-7;
end
pam_sig = pam_sig/sqrt(21);
@@ -233,13 +257,19 @@ classdef PAMmapper
case 2
% 2-ASK
data_out=comp_real(:,:,1);
if ~obj.eth_style
data_out=comp_real(:,:,1);
else
data_out=abs(comp_real(:,:,1)-1);
end
case 4
% 4-ASK
data_out=[comp_real(:,:,2); ones(s1,s2) - comp_real(:,:,1) + comp_real(:,:,3)];
if ~obj.eth_style
data_out=[comp_real(:,:,2); ones(s1,s2) - comp_real(:,:,1) + comp_real(:,:,3)];
else
data_out= [(data_in>=0); (abs(data_in)<=1)];
end
case 6
@@ -262,13 +292,18 @@ classdef PAMmapper
data_out(m:m+4) = bitget(dd_idx-1,5:-1:1);
m = m+5;
end
case 8
% 8-ASK
if ~obj.eth_style
data_out=[comp_real(:,:,4);
comp_real(:,:,1)-comp_real(:,:,3)+comp_real(:,:,5)-comp_real(:,:,7);
1-comp_real(:,:,2)+comp_real(:,:,6)];
else
data_out = [(data_in>=0); (abs(data_in)>(4/sqrt(21))); (abs(data_in)>=(2/sqrt(21)))&(abs(data_in)<=(6/sqrt(21)))];
end
case 16
% 16-ASK

27
Classes/01_transmit/PAMsource.m
View File

@@ -83,6 +83,8 @@ classdef PAMsource
%%%%% MOVE-IT PRMS %%%%
state = struct();
para = struct();
@@ -103,6 +105,8 @@ classdef PAMsource
para.reset_prms = 0;
para.method = 1;
data_in = [];
global loop;
loop = 0;
@@ -151,12 +155,16 @@ classdef PAMsource
sym_max = max(symbols.signal);
end
% symbols.move_it_spectrum("fignum",222,"displayname","Symbols only");
% symbols.spectrum("fignum",222,"displayname","Symbols only","normalizeTo0dB",1);
% symbols.spectrum("fignum",111,"displayname","1) RAW SIGNAL");
pulsf = Moveit_wrapper('pulsef');
pulsf.para.alpharacos = obj.pulseformer.alpha;
pulsf.para.f_sym = obj.fsym;
pulsf.para.fs = obj.fsym;
pulsf.para.pulse = 'racos';
pulsf.process(symbols);
%%%%% Pulse-forming %%%%%%
if obj.applypulseform
digi_sig = obj.pulseformer.process(symbols);
@@ -164,10 +172,13 @@ classdef PAMsource
digi_sig = symbols;
end
%%%%% Re-sample to f DAC %%%%%%
digi_sig = digi_sig.resample("fs_in",digi_sig.fs,"fs_out",obj.fs_out,"n",10,"beta",5);
% digi_sig.spectrum("fignum",112,"displayname","2) after pulseforming");
%%%%% Re-sample to f DAC %%%%%%
n = 10;
digi_sig = digi_sig.resample("fs_in",digi_sig.fs,"fs_out",obj.fs_out,"n",n,"beta",5);
% digi_sig.spectrum("fignum",111,"displayname",['3) Shaped + Resampled; n: ',num2str(n)]);
% digi_sig.spectrum("fignum",111,"displayname","after pulseforming");
%%%%% Hard clip digital signal to PAM range before DAC %%%%%%
if obj.applyclipping

112
Classes/01_transmit/Pulseformer.m
View File

@@ -3,13 +3,15 @@ classdef Pulseformer
% Detailed explanation goes here
properties(Access=public)
fdac
end
properties(Access=private)
properties(Access=public)
fdac
fsym
pulse
pulselength
rrcalpha
alpha
matched
end
methods (Access=public)
@@ -20,9 +22,10 @@ classdef Pulseformer
arguments
options.fdac double
options.fsym double
options.pulse pulseform = pulseform.rrc
options.pulse pulseform = pulseform.rc
options.pulselength double {mustBeInteger} = 32
options.rrcalpha double = 0.05
options.alpha double = 0.05
options.matched = 0;
end
%
@@ -86,49 +89,72 @@ classdef Pulseformer
sps= up;
end
if obj.pulse == pulseform.rrc
%Bau das Filter (hier rrc)
racos_len = obj.pulselength*2;
alpha = obj.rrcalpha;
h = rcosdesign(alpha,racos_len,sps,"normal");
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 = 1./h;
end
manual_cyclic_convolution = 1;
upfirdn_convolution = 0;
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:dn:end);
end
end
if upfirdn_convolution
%Apply Filter using Matlab build in fctn.
h = rcosdesign(obj.alpha,racos_len,sps,"normal");
h = h./ max(h);
data_out_ = upfirdn(data_in,h,up,dn);
%cut signal, which is longer due to fir filter
st = round(up/dn*racos_len/2)+1; %we need to cut y_out
en = round(st + (length(data_in)*up/dn));
data_out_ = data_out_(st:en);
end
% Apply filter the long way (from move_it)
% block length in samples
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
% cyclic convolution
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:dn: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
% %Apply Filter using Matlab build in fctn.
% h = rcosdesign(alpha,racos_len,sps);
% h = h./ max(h);
%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);
%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;

2
Classes/Warehouse_class/classes/DataStorage.m
View File

@@ -89,7 +89,7 @@ classdef DataStorage < handle
for p = 1:numel(obj.fn)
name = obj.fn(p);
values = obj.inputParams.(name);
obj.parameter.(name) = Parameter(name,values);
obj.parameter.(name) = StorageParameter(name,values);
end
end

44
Classes/Warehouse_class/classes/Parameter.m
View File

@@ -1,44 +0,0 @@
classdef Parameter < handle
%PARAMETER Summary of this class goes here
% Detailed explanation goes here
properties
name
values
indices
length
getPhysForIndex = dictionary;
getIndexForPhys = dictionary;
end
methods
function obj = Parameter(name, values)
%PARAMETER Construct
obj.name = name;
obj.values = values;
obj.indices = [1:numel(obj.values)];
obj.length = length(obj.values);
obj = obj.buildIndexForPhysDict();
obj = obj.buildPhysForIndexDict();
end
function obj = buildPhysForIndexDict(obj)
%take Index
%return Phys
for idx = 1:numel(obj.values)
obj.getPhysForIndex(idx) = obj.values(idx);
end
end
function obj = buildIndexForPhysDict(obj)
%take Phys
%return Index
for idx = 1:numel(obj.values)
obj.getIndexForPhys(obj.values(idx)) = idx;
end
end
end
end

3
Datatypes/pulseform.m
View File

@@ -1,7 +1,8 @@
classdef pulseform < int32
enumeration
rrc (1)
rc (1) %raised cosine (use this at Tx without matched filter)
rrc (2) %root raised cosine (usually with matched filter)
end
end

13
Functions/EQ_structures/vnle.m
View File

@@ -22,6 +22,7 @@ function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options)
tx_bits
options.precode_mode db_mode
options.showAnalysis = 0
options.eth_style = 0;
end
%FFE or VNLE
@@ -46,7 +47,7 @@ function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options)
tx_symbols_precoded = Duobinary().encode(tx_symbols);
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
tx_bits = PAMmapper(M,0).demap(tx_symbols_precoded);
tx_bits = PAMmapper(M,0,"eth_style",options.eth_style).demap(tx_symbols_precoded);
case db_mode.db_discard
@@ -64,23 +65,31 @@ function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options)
end
rx_bits = PAMmapper(M,0).demap(eq_signal_hd);
rx_bits = PAMmapper(M,0,"eth_style",options.eth_style).demap(eq_signal_hd);
[~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
[evm_total,evm_lvl] = calc_evm(eq_signal_sd,tx_symbols);
[inf_rate] = calc_air(eq_signal_sd,tx_symbols,"skip_front",10000,"skip_end",10000);
snr(eq_signal_sd.signal,eq_noise.signal)
eq_package.ber_vnle = ber;
eq_package.evm_total = evm_total;
eq_package.evm_lvl = evm_lvl;
eq_package.inf_rate_vnle = inf_rate;
if options.showAnalysis
fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl);
fprintf('VNLE BER: %.2e \n',ber);
showEQNoiseSNR(tx_symbols,eq_noise,"displayname",'SNR after VNLE','fignum',301);
showLevelHistogram(eq_signal_sd,tx_symbols,"fignum",302);
end
end

2
Functions/EQ_structures/vnle_postfilter_mlse.m
View File

@@ -99,7 +99,7 @@ end
fprintf('MLSE BER: %.2e \n',ber_mlse);
showEQNoisePSD(eq_noise,"fignum",336,"displayname",'VNLE+DFE','postfilter_taps',pf_.coefficients);
showEQNoisePSD(eq_noise,"fignum",336,"displayname",'Residual Noise after VNLE','postfilter_taps',pf_.coefficients);
rx_signal.spectrum("normalizeTo0dB",1,"fignum",337,"displayname",'Rx Signal');

25
Functions/EQ_visuals/showEQNoiseSNR.m
View File

@@ -1,7 +1,7 @@
function showEQNoiseSNR(tx_signal, rx_signal, options)
function showEQNoiseSNR(eq_signal, noise_signal, options)
arguments
tx_signal
rx_signal
eq_signal
noise_signal
options.fs_tx
options.fs_rx
options.fignum (1,1) double = NaN % Default to NaN if not provided
@@ -16,13 +16,13 @@ end
fig = figure(options.fignum); % Use the specified figure number
end
if isa(tx_signal,'Signal')
options.fs_tx = tx_signal.fs;
tx_signal = tx_signal.signal;
if isa(eq_signal,'Signal')
options.fs_tx = eq_signal.fs;
eq_signal = eq_signal.signal;
end
if isa(rx_signal,'Signal')
options.fs_rx = rx_signal.fs;
rx_signal = rx_signal.signal;
if isa(noise_signal,'Signal')
options.fs_rx = noise_signal.fs;
noise_signal = noise_signal.signal;
end
@@ -37,10 +37,10 @@ end
% Ensure the figure is ready before calling spectrum
title('SNR of received Signal')
fft_length = 2^(nextpow2(length(tx_signal))-7);
fft_length = 2^(nextpow2(length(eq_signal))-7);
[s_lin,w] = pwelch(tx_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_tx,"centered","psd","mean");
[n_lin,w] = pwelch(rx_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_rx,"centered","psd","mean");
[s_lin,w] = pwelch(eq_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_tx,"centered","psd","mean");
[n_lin,w] = pwelch(noise_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_rx,"centered","psd","mean");
w = w.*1e-9;
snr_dbm = 10*log10(s_lin./n_lin);
@@ -48,6 +48,7 @@ end
% figure(231)
hold on
plot(w,snr_dbm,'DisplayName','SNR','LineWidth',0.5,'Color',options.color);
yline(mean(snr_dbm),'HandleVisibility','off','Color',options.color);
xlabel("Frequency in GHz");
edgetick = 2^(nextpow2(options.fs_tx*1e-9));

BIN
projects/.DS_Store vendored
View File

Binary file not shown.

8
projects/IMDD_base_system/imdd_it.m
View File

@@ -3,13 +3,13 @@
% db = DBHandler("pathToDB",[basePath,'silas_labor.db']);
if 1
uloops = struct;
uloops.precomp = [0,1];
uloops.db_precode = [0,1];
uloops.bitrate = [300,330,360,390,420,450,480].*1e9; %[300,330,360,390,420,450,480]
uloops.precomp = [0];
uloops.db_precode = [0];
uloops.bitrate = [330].*1e9; %[300,330,360,390,420,450,480]
% uloops.bitrate = 390e9;
% uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4];
uloops.laser_wavelength = [1293];
uloops.M = [4,6,8];
uloops.M = [4];
uloops.link_length = [2]; % 1,2,3,5,6,8,10
wh = DataStorage(uloops);
wh.addStorage("ber");

8
projects/IMDD_base_system/imdd_model.m
View File

@@ -135,7 +135,7 @@ if simulation_mode
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rcalpha);
[Digi_sig,Symbols,Tx_bits] = PAMsource(...
"fsym",fsym,"M",M,"order",17,"useprbs",1,...
"fsym",fsym,"M",M,"order",16,"useprbs",0,...
"fs_out",fdac,...
"applyclipping",0,"clipfactor",1.5,...
"applypulseform",apply_pulsef,"pulseformer",Pform,...
@@ -268,7 +268,7 @@ for occ = 1:proc_occ
Scpe_sig = Scpe_sig.resample("fs_out",2*fsym);
%%%%%% Sync Rx signal with reference %%%%%%
% [Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym);
[Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym);
Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig);
@@ -284,10 +284,10 @@ for occ = 1:proc_occ
mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
% %%%%% VNLE + DFE %%%%
if 0
if 1
eq_vnle_dfe = EQ("Ne",vnle_order,"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",1,"ideal_dfe",0);
[result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",0);
[result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",1);
vnle_dfe_package{occ} = result;