diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 48ddbf0..fa4c44a 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -237,12 +237,16 @@ classdef Signal arguments obj Signal - options.fs_in double + options.fs_in double = obj.fs options.fs_out double options.n double = 10; options.beta double = 5; end + if options.fs_in ~= obj.fs + warning('The signals fs is different from the given fs_in while it should be the same.'); + end + obj.signal = resample(obj.signal,options.fs_out,options.fs_in,options.n,options.beta); desc = ['resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ]; @@ -265,7 +269,7 @@ classdef Signal % spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname); N = 2^(nextpow2(length(obj.signal))-8); - + [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" @@ -279,7 +283,7 @@ classdef Signal edgetick = 2^(nextpow2(obj.fs*1e-9)); % xticks([-edgetick:16:edgetick]); xlim([-244, 244]) - ylim([-120,-0]); + ylim([-120,10]); yticks([-200:10:10]); legend @@ -301,7 +305,7 @@ classdef Signal pow = pow / 50; end pow = 10*log10(pow)+30; %dbm - + case power_notation.mW pow = pow .* 1e3; %mW @@ -378,14 +382,23 @@ classdef Signal end %% - function [obj,delay_n] = delay(obj,options) + function [obj] = delay(obj,delay,options) arguments obj Signal - options.delay_samples = 0 + delay double = 0 + options.mode delay_mode = delay_mode.samples end - obj.signal=delayseq(obj.signal,options.delay_samples); + if options.mode == delay_mode.samples + + obj.signal=delayseq(obj.signal,delay); + + elseif options.mode == delay_mode.time + + obj.signal=delayseq(obj.signal,delay,obj.fs); + + end end @@ -415,7 +428,7 @@ classdef Signal end % delay by lagging samples - obj = obj.delay("delay_samples",-D); + obj = obj.delay(-D,'mode','samples'); cuts = obj.length-(options.reference.length*q); @@ -491,7 +504,7 @@ classdef Signal ppeak(i) = maxA - (difference/histpoints*loc(i)); end - + if isa(obj,'Opticalsignal') er=10*log10(ppeak(1)/ppeak(end)); @@ -633,31 +646,7 @@ classdef Signal 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] @@ -687,24 +676,53 @@ classdef Signal '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 + try + 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'); - % Create third annotation box for Vmax - if ~oma - thirdboxstring = ['ER (db):',num2str(er),' dB']; - else - thirdboxstring = ['OMA outer:',num2str(er),' V']; - end + [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',':'); - annotation('textbox', boxPosition, ... - 'String',thirdboxstring , ... - 'BackgroundColor', boxColor, ... - 'EdgeColor', boxEdgeColor, ... - 'LineStyle', boxLineStyle, ... - 'FontWeight', boxFontWeight, ... - 'HorizontalAlignment', 'center'); + 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 + + % 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'); + end yticks(linspace(0,histpoints,16)); y_tickstring = sprintfc('%.2f', y_tickstring); diff --git a/Classes/01_transmit/ChannelFreqResp.m b/Classes/01_transmit/ChannelFreqResp.m new file mode 100644 index 0000000..d927257 --- /dev/null +++ b/Classes/01_transmit/ChannelFreqResp.m @@ -0,0 +1,406 @@ +classdef ChannelFreqResp < handle + % Linear pre-compensation of Channel effects. This module has two modes: + + % A) In acquire mode it sends a real OFDM singal with all subcarriers assigned in + % order to acquire later the system's frequency response after the + % signal passed the system... + % ChannelFreqResp.SendOFDM + + % B) In the "not acquire" mode it loads the frequency response generated by the freq_res module, + % calculates the inverse frequency response with the right length matched + % with the signal length and distorts the signal before tranmission + % give channel indices for measure mode + + %TYPICAL FLOW: + + % referencesignal = ChannelFreqResp.buildOFDM() + + % referencesignal -> SYSTEM -> measuredsignal + + % ChannelFreqResp.estimate(measuredsignal) + + properties(Access=public) + Nacq + Navg + Ncp + + f_ref %fs of transmitted dmt sig + f_observe %fs of the observed signal + + symlen + seqlen + + refsig + H + H_all + H_apply + + Nfft + df + faxis + end + + methods (Access=public) + + function obj = ChannelFreqResp(options) + %NAME Construct an instance of this class + % Detailed explanation goes here + + arguments + options.Nacq = 4096; + options.Navg = 30; + options.Ncp = 100; + options.f_ref; + options.f_observe; + end + + % + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + + obj.Nfft = 2*obj.Nacq + 1 ; + obj.df = obj.f_ref/obj.Nfft ; % calculate frequency grid + obj.faxis = (0:obj.Nfft-1)*obj.df ; + + end + + function output = buildOFDM(obj) + % generate OFDM or DMT signal that is send over the channel + + obj.symlen = (obj.Nacq*2+obj.Ncp+1); + + obj.seqlen = obj.symlen*obj.Navg; + + obj.refsig = obj.genRNDDMT(15)'; + + obj.refsig = Informationsignal(obj.refsig,"fs",obj.f_ref); + + output = obj.refsig; + + end + + function output = estimate(obj, data_in, options) + arguments + obj + data_in + options.save logical = false + options.savePath = ""; + options.fileName = ""; + end + % use the transmitted DMT and the received DMT to get the transfer function + + obj.Nfft = 2*obj.Nacq + 1 ; + + %resample from fadc to fdac (fdac => fs of reference signal) + data_in = data_in.resample("fs_out",obj.f_ref); + + %est from dmt + obj.H_all = obj.estHfromDMT(data_in.signal,obj.refsig.signal); + + obj.H = mean(obj.H_all,1); + + output = obj.H; + + if options.save + obj.save('fileName',options.fileName,'savePath',options.savePath); + end + + end + + function Target = precomp(obj,Target,options) + % apply the freq response + arguments + obj + Target + options.fileName = '' + options.loadPath = '' + options.maxampdb + end + + if isempty(obj.H) + obj.load("fileName",options.fileName,"loadPath",options.loadPath); + end + + H_inv = 1./obj.H; + + nH = find(~isnan(H_inv),1,'last'); %last value that is not nan + H_inv(nH:end)=H_inv(nH); %replace everything after first NaN with the last non-NaN value + + fstarget = Target.fs; + + % Build new frequencie axis (with current fs) + fnew = linspace(0,fstarget/2,length(Target)/2+1); + fnew = fnew(2:end-1); + + % Old frequency axis (should be much coarser) + idx_old = find((obj.faxis > 0) .* (obj.faxis < fstarget/2)); %positions of all Frequencies smaller than fs/2 + int_fold = obj.faxis(idx_old); %old frequencies from 0 to fs/2 + + idx_new = find(fnew <= int_fold(end)); %positions of all Frequencies smaller than fs/2 + int_fnew = fnew(idx_new);%new frequencies from 0 to fs/2 + + % interpolate the frequency response that had a coarse frequency resolution (e.g. 256 bins) to the current frequency resolution (e.g. 21843 bins) + iH = interp1(int_fold, real(H_inv(idx_old)) ,fnew, 'linear') + 1i*interp1(int_fold, imag(H_inv(idx_old)) ,fnew, 'linear'); + + % set all NaN values to the fist/ last non-NaN value + nH = find(~isnan(iH),1,'first'); + iH(1:nH)=iH(nH); + nH = find(~isnan(iH),1,'last'); + iH(nH:end) =iH(nH); + + %smoothing takes time and sometimes the result looks odd, + %however the performance is most of the time better + smoothing = 1; + if smoothing + iH = smooth(fnew,iH,0.1,'loess')'; + end + + % multiply the complex frequency responce with the phase + % angle (-pi,pi) at lowest frequency + iH = iH.*exp(-1j*angle(iH(1))); % to be checked (<- not from silas, so what needs to be checked?) + + % Phase difference between lowest and hiughest frequency + % component -> but what is this for? dPhase is not used... + noncausal = 0; + if noncausal + dPhase = angle(iH(end))-angle(iH(1)); + end + + % normalize complex freq. resp. by magnitude at the first + % five frequencies -> should be the vaue at f=0=DC component? + iH = iH./mean(abs(iH(1:100))); %why 1:5?? + + % set maximum amplification + % set als values higher than hmax to hmax and keep the + % phase information by multiplication with respective + % corresponding phase angles + maxamp_lin = 10^(options.maxampdb/20); + + iH(abs(iH)>maxamp_lin) = maxamp_lin.*exp(1j*angle(iH(abs(iH)>maxamp_lin))); + + % it could be helpful to clip at linear 1 (=to keep comp from attenuating the signal) + % iH(abs(iH)<1) = 1.*exp(1j*angle(iH(abs(iH)<1))); + + if 0 + figure(7);hold on;plot(fnew,20*log10(abs(iH))) + end + H_inv = [iH(1) iH fliplr(conj(iH)) conj(iH(1))]; + + obj.H_apply = H_inv; + + Target.signal = real((ifft( ( fft(real( Target.signal )) .* H_inv' ) ))); + + + end + + function plot(obj) + + figure(55551); + clf; + + Havg = obj.H; + + %1) + subplot(4,1,1);hold all;box on;title('Magnitude Freq. Response'); + plot(obj.faxis/1e9, 20*log10(abs(obj.H_all)),'linewidth',0.1,'LineStyle','-','Color','#808080') ; + xlim([0.2 .5*max(obj.faxis)*1e-9]); + plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2); + grid on; + + %2) + subplot(4,1,2); hold all; box on; title('Phase Freq. Response'); + plot(obj.faxis/1e9, (angle(obj.H_all)),'linewidth',0.1,'LineStyle','-','Color','#808080') ; + plot(obj.faxis/1e9, (angle(Havg)),'LineWidth',2) ; + xlim([0.2 .5*max(obj.faxis)*1e-9]); + grid on; + + %normalize / remove attenuation + Havg = Havg./mean(Havg(2:10)); + + %3) + subplot(4,1,3); hold all; box on; title('Inverse Magnitude Freq. Response'); + plot(obj.faxis/1e9, 20*log10(abs(1./Havg)),"LineWidth",2,"Color",[0.3467 0.5360 0.6907]) ; + xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on; + ylim([-1 15]); + hold on; + yline(3,'LineWidth',2,'LineStyle','--'); + + %4) + subplot(4,1,4); hold all; box on; title('Inverse Phase Freq. Response'); + plot(obj.faxis/1e9, (angle(1./Havg)),"LineWidth",2,"Color",[0.3467 0.5360 0.6907]) ; + xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on; + + + end + + function save(obj,options) + arguments + obj + options.fileName = '' + options.savePath = '' + end + + % Check if the fileName was provided + if isempty(char(options.fileName)) + % If no file name is provided, prompt the user to enter a file name + [options.fileName, filePath] = uiputfile('*.mat', 'Save As'); + + % If the user cancels the dialog, fileName and filePath will be 0 + if isequal(options.fileName, 0) || isequal(filePath, 0) + disp('Save operation cancelled.'); + return; + end + + % Update the savePath with the directory chosen by the user + options.savePath = filePath; + end + + % Check if the savePath was provided + if isempty(char(options.savePath)) + % If no path is provided, open a UI window to select the path + options.savePath = uigetdir('', 'Select a folder to save the file'); + + % If the user cancels the dialog, savePath will be 0 + if options.savePath == 0 + disp('Save operation cancelled.'); + return; + end + else + % If a path is provided, validate it + if ~isfolder(options.savePath) + error('The specified path does not exist.'); + end + end + + % Construct the full file path + fullFileName = fullfile(options.savePath, options.fileName); + + % Save the data to the specified file + save(fullFileName, 'obj'); + fprintf('Frequency response information successfully saved to %s\n', fullFileName); + + end + + function data = load(obj, options) + % Function to load data from a specified file and path. + arguments + obj + options.fileName = '' + options.loadPath = '' + end + + % Check if the fileName was provided + if isempty(char(options.fileName)) + % If no file name is provided, open a UI window to select the file + [options.fileName, options.loadPath] = uigetfile('*.mat', 'Select a file to load'); + + % If the user cancels the dialog, fileName and loadPath will be 0 + if isequal(options.fileName, 0) || isequal(options.loadPath, 0) + disp('Load operation cancelled.'); + data = []; + return; + end + end + + % If loadPath is not provided or empty, use the current folder + if isempty(char(options.loadPath)) + options.loadPath = pwd; + else + % Validate the path if provided + if ~isfolder(options.loadPath) + error('The specified path does not exist.'); + end + end + + % Construct the full file path + fullFileName = fullfile(options.loadPath, options.fileName); + + % Check if the file exists + if ~isfile(fullFileName) + fullFileName = fullfile(options.loadPath, [char(options.fileName),'.mat']); + if ~isfile(fullFileName) + error('The specified file does not exist.'); + end + end + + % Load the data from the specified file + loadedData = load(fullFileName); + + % Replace whole obj here.. is this save or unsave?! + fn = fieldnames(loadedData.obj); + for n = 1:numel(fn) + try + obj.(fn{n}) = loadedData.obj.(fn{n}); + end + end + + fprintf('Frequency response information successfully loaded from %s\n', fullFileName); + end + + end + + methods (Access=private) + % Cant be seen from outside! So put all your functions here that can/ + % shall not be called from outside + + function rOFDM = genRNDDMT(obj,randkey) + + rOFDM = NaN(1,(2*obj.Nacq+obj.Ncp+1)*obj.Navg); + + s = RandStream('mt19937ar','Seed',randkey,'NormalTransform','Polar'); + + ref = 2 * round(rand(s,obj.Nacq, obj.Navg)) - 1 ; % Navg Random BPSK sequences + + ofdm = [ones(1, obj.Navg) ; ref ; conj(ref(end:-1:1,:)) ] ; % DMT + + ofdm = ifft(ofdm) ; % Navg real OFDM sequences + + ofdm = [ofdm(end-obj.Ncp+1 : end, :) ; ofdm] ; % Add cyclic prefix + + rOFDM = reshape(ofdm, 1, size(ofdm,1)*size(ofdm,2)) ; + + rOFDM = rOFDM./max(abs(rOFDM(:))); + + end + + function [rH] = estHfromDMT(obj, data_in, ref_in) + + %! Dont (circ)shift the signal here as this would remove the phase information! + %tested with a butterworth filter this exactly reconstructs the + %phase and the magnitude. However, the option is here + estimatephase = 1; + + if ~estimatephase + % 0. cross-correlate the received signal with its reference to extract the periods + corr = abs(ifft( fft(data_in(1:length(ref_in))).* conj(fft(ref_in)) )) ; + + % find max + [~, peak] = max(corr) ; + peak=max(1,peak-1); + + data_in = circshift(data_in,-peak) ; + %Y = data_in(peak:peak+(Nfft+obj.Ncp)*obj.Navg-1) ; + end + + % 1. Reshape signal to a matrix to support noise averaging + Nfft = 2*obj.Nacq + 1 ; + + Y = reshape(data_in, Nfft+obj.Ncp, obj.Navg).' ; + + X = reshape(ref_in, Nfft + obj.Ncp, obj.Navg).' ; + + % 2. Remove cyclic prefix and apply FFT transformation + Y = fft(Y(:, obj.Ncp+1 : obj.Ncp + Nfft), [], 2) ; + X = fft(X(:, obj.Ncp+1 : obj.Ncp + Nfft), [], 2) ; + + % 4. Caclulate the frequency response using H = Y/X + rH = Y./X ; + + end + + end +end diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index ed015fe..2a66d6d 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -24,7 +24,9 @@ classdef PAMmapper if isa(signal_in,'Signal') signal_in.signal = obj.map_(signal_in.signal); - signal_in = signal_in.normalize("mode","rms"); + + % signal_in = signal_in.normalize("mode","rms"); + signal_in = signal_in.logbookentry(); out = signal_in; else @@ -42,7 +44,7 @@ classdef PAMmapper function pam_sig = map_(obj,bitpattern) switch obj.M - case 1 + case 2 % 2-ASK: BPSK / OOK pam_sig=bitpattern(:,1); @@ -87,6 +89,7 @@ classdef PAMmapper end pam_sig = pam_sig/sqrt(21); + case 16 % 16-ASK: x1 = bitpattern(:,1); diff --git a/Classes/01_transmit/PAMsource.m b/Classes/01_transmit/PAMsource.m index 8df05af..e72033c 100644 --- a/Classes/01_transmit/PAMsource.m +++ b/Classes/01_transmit/PAMsource.m @@ -9,6 +9,8 @@ classdef PAMsource fsym randkey + db_precode + mrds_code mrds_blocklength @@ -33,6 +35,8 @@ classdef PAMsource options.fsym = 112e9; options.randkey = 0; + options.db_precode = 0; + options.mrds_code = 0; options.mrds_blocklength = 512; @@ -90,6 +94,11 @@ classdef PAMsource symbols = PAMmapper(obj.M,0).map(bits); symbols.fs = obj.fsym; + + if obj.db_precode + symbols = Duobinary().precode(symbols); + symbols = Duobinary().encode(symbols); + end if obj.mrds_code symbols = MRDS_coding("blocklength",obj.mrds_blocklength).encode(symbols); diff --git a/Classes/04_DSP/A1_scheme.m b/Classes/04_DSP/A1_scheme.m deleted file mode 100644 index 0a76c1c..0000000 --- a/Classes/04_DSP/A1_scheme.m +++ /dev/null @@ -1,23 +0,0 @@ -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/Coding/Duobinary.m b/Classes/04_DSP/Coding/Duobinary.m new file mode 100644 index 0000000..8049078 --- /dev/null +++ b/Classes/04_DSP/Coding/Duobinary.m @@ -0,0 +1,245 @@ +classdef Duobinary + %Duobinary Coding + % should work + + properties(Access=public) + + + end + + methods (Access=public) + + function obj = Duobinary() + %NAME Construct an instance of this class + % Detailed explanation goes here + + arguments + + end + + % % + % fn = fieldnames(options); + % for n = 1:numel(fn) + % try + % obj.(fn{n}) = options.(fn{n}); + % end + % end + + % do more stuff + + end + + function signalclass = precode(~,signalclass) + + data = signalclass.signal; + u = unique(data); + M = numel(u); + + %make unipolar + if M == 4 + data = data .* sqrt(5); + elseif M == 6 + data = data .* sqrt(10); + elseif M == 8 + data = data .* sqrt(21); + elseif M == 16 + data = data .* sqrt(85); + warning('Check if PAM16 implementation, mapping and scaling is correct!') + end + + data = round(data); + b = min(data); + data = data - b; + data = data ./ 2; + + assert(isequal((0:M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar + + % Pre coding + bk = zeros(numel(data),1); + + for k = 1:numel(data)-1 + bk(k+1) =mod(data(k)-bk(k),M); + end + + %make bipolar + bk = bk .* 2; + bk = bk + b; + + if M == 4 + bk = bk ./ sqrt(5); + elseif M == 6 + bk = bk ./ sqrt(10); + elseif M == 8 + bk = bk ./ sqrt(21); + end + signalclass.signal = bk; + + assert(isequal(unique(signalclass.signal),u),'Check Duobinary Precoding'); %seems the signal is not the same as before + + end + + function signalclass = encode(~,signalclass) + + data = signalclass.signal; + u = unique(data); + M = numel(u); + + %make unipolar + if M == 4 + data = data .* sqrt(5); + elseif M == 6 + data = data .* sqrt(10); + elseif M == 8 + data = data .* sqrt(21); + elseif M == 16 + data = data .* sqrt(85); + warning('Check if PAM16 implementation, mapping and scaling is correct!') + end + + data = round(data); + b = min(data); + data = data - b; + data = data ./ 2; + + assert(isequal((0:M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar + + % duobinary coding (1+D) + coeff = [1,1]; + + data = conv(data,coeff,"same"); + + %make bipolar + data = (data-round(mean(data),1)); + + if M == 4 + data = data ./ sqrt(2.5); % 7-level constellation weighted with probability after DB code i.e. mean([-3 3 -2 -2 2 2 -1 -1 -1 1 1 1 0 0 0 0].^2) = 2.5 --> sqrt(2.5) == rms(constellation) + elseif M == 6 + data = data ./ sqrt(5.8); + elseif M == 8 + data = data ./ sqrt(10.5); % 15-level constellation weighted with probability after DB code i.e. + end + + signalclass.signal = data; + + % Code to evaluate the s in sqrt(s): + % M=6; + % c = (0:(2*M)-2) - ((2*M)-2)/2; + % s = 0; + % for i = 0:(2*M)-2 + % p(i+1) = M-abs(i-(M-1)); + % s = s+p(i+1)*c(i+1).^2; + % end + % s = s/M^2; + + end + + function signalclass = decode(~,signalclass) + + data = signalclass.signal; + u = unique(data); + I = numel(u); %number of duobinary coded const. points + M = (I+1)/2; %PAM-M order + + %make unipolar + if I == 7 + data = data .* sqrt(2.5); + elseif I == 11 + %todo + data = data .* sqrt(5.8); + warning('Check if PAM16 implementation, mapping and scaling is correct!') + elseif I == 15 + data = data .* sqrt(10.5); + elseif I == 16 + warning('Check if PAM16 implementation, mapping and scaling is correct!') + end + + data = round(data); + b = min(data); + data = data - b; + + data = round(data); + + data = mod(data,M); + + %make bipolar + data = data .* 2; + data = data - round(mean(data)); + + if M == 4 + data = data ./ sqrt(5); + elseif M == 6 + data = data ./ sqrt(10); + elseif M == 8 + data = data ./ sqrt(21); + end + signalclass.signal = data; + + end + + function signalclass = feedbackdetector(~,signalclass) + data = signalclass.signal; + u = unique(data); + I = numel(u); %number of duobinary coded const. points + M = (I+1)/2; %PAM-M order + + %make unipolar + if I == 7 + data = data .* sqrt(2.5); + elseif I == 11 + %todo + data = data .* sqrt(5.8); + warning('Check if PAM16 implementation, mapping and scaling is correct!') + elseif I == 15 + data = data .* sqrt(10.5); + elseif I == 16 + warning('Check if PAM16 implementation, mapping and scaling is correct!') + end + + data = round(data); + b = min(data); + data = data - b; + + data = round(data); + + data_out = zeros(length(data),1); + const = 0:M-1; + %%FALSCH! + for k = 1:length(data)-1 + + d_ = data(k+1) - data_out(k); + + [~,b] = min(abs(d_-[0:M-1])); + + data_out(k+1) = const(b); + + end + + %make bipolar + data_out = data_out .* 2; + data_out = data_out - round(mean(data_out)); + + if M == 4 + data_out = data_out ./ sqrt(5); + elseif M == 6 + data_out = data_out ./ sqrt(10); + elseif M == 8 + data_out = data_out ./ sqrt(21); + end + + signalclass.signal = data_out; + + + + + 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/04_DSP/MRDS_coding.m b/Classes/04_DSP/Coding/MRDS_coding.m similarity index 100% rename from Classes/04_DSP/MRDS_coding.m rename to Classes/04_DSP/Coding/MRDS_coding.m diff --git a/Classes/04_DSP/EQ.m b/Classes/04_DSP/Equalizer/EQ.m similarity index 100% rename from Classes/04_DSP/EQ.m rename to Classes/04_DSP/Equalizer/EQ.m diff --git a/Classes/04_DSP/EQ_copy.m b/Classes/04_DSP/Equalizer/EQ_copy.m similarity index 100% rename from Classes/04_DSP/EQ_copy.m rename to Classes/04_DSP/Equalizer/EQ_copy.m diff --git a/Classes/04_DSP/EQ_silas.m b/Classes/04_DSP/Equalizer/EQ_silas.m similarity index 97% rename from Classes/04_DSP/EQ_silas.m rename to Classes/04_DSP/Equalizer/EQ_silas.m index 42470af..a63db84 100644 --- a/Classes/04_DSP/EQ_silas.m +++ b/Classes/04_DSP/Equalizer/EQ_silas.m @@ -1,486 +1,486 @@ -classdef EQ_silas < 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 - d_out %decision output - - % 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(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,symbols_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'; - - - %change sampling frequency of outgoing signal - signalclass_in.fs = reference_signalclass_in.fs; - - % append to logbook - lbdesc = ['EQ von Silas ist gelaufen ']; - signalclass_in = signalclass_in.logbookentry(lbdesc); - - symbols_out = signalclass_in; - symbols_out.signal = obj.d_out; - - % 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); - - 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 - - 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 = 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; - - 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).'; - - %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) = x_d.'* coeff; - - %Decision 1 - [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point - d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx); - - y_1(m,symbol_idx) = y(m); % after 1st iteration - - %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(m) * conj(x_d))); - - % 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(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 - - - - %% - obj.y_out = (circshift( y.' ,-(obj.delay))).'; - obj.d_out = d_hat(1:2:end); - - % 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 - - %% 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 - - 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 - - - - - - - - - - - - - - - - - - - - - - +classdef EQ_silas < 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 + d_out %decision output + + % 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(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,symbols_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'; + + + %change sampling frequency of outgoing signal + signalclass_in.fs = reference_signalclass_in.fs; + + % append to logbook + lbdesc = ['EQ von Silas ist gelaufen ']; + signalclass_in = signalclass_in.logbookentry(lbdesc); + + symbols_out = signalclass_in; + symbols_out.signal = obj.d_out; + + % 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); + + 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 + + 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 = 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; + + 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).'; + + %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) = x_d.'* coeff; + + %Decision 1 + [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point + d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx); + + y_1(m,symbol_idx) = y(m); % after 1st iteration + + %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(m) * conj(x_d))); + + % 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(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 + + + + %% + obj.y_out = (circshift( y.' ,-(obj.delay))).'; + obj.d_out = d_hat(1:2:end); + + % 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 + + %% 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 + + 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_ofc.m b/Classes/04_DSP/Equalizer/EQ_silas_ofc.m similarity index 100% rename from Classes/04_DSP/EQ_silas_ofc.m rename to Classes/04_DSP/Equalizer/EQ_silas_ofc.m diff --git a/Classes/04_DSP/EQ_silas_plain.m b/Classes/04_DSP/Equalizer/EQ_silas_plain.m similarity index 100% rename from Classes/04_DSP/EQ_silas_plain.m rename to Classes/04_DSP/Equalizer/EQ_silas_plain.m diff --git a/Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m b/Classes/04_DSP/Equalizer/EQ_silas_sliding_window_dc_removal.m similarity index 100% rename from Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m rename to Classes/04_DSP/Equalizer/EQ_silas_sliding_window_dc_removal.m diff --git a/Classes/04_DSP/FFE.m b/Classes/04_DSP/Equalizer/FFE.m similarity index 100% rename from Classes/04_DSP/FFE.m rename to Classes/04_DSP/Equalizer/FFE.m diff --git a/Classes/04_DSP/FFE_DCremoval.m b/Classes/04_DSP/Equalizer/FFE_DCremoval.m similarity index 100% rename from Classes/04_DSP/FFE_DCremoval.m rename to Classes/04_DSP/Equalizer/FFE_DCremoval.m diff --git a/Classes/04_DSP/FFE_DFE.m b/Classes/04_DSP/Equalizer/FFE_DFE.m similarity index 100% rename from Classes/04_DSP/FFE_DFE.m rename to Classes/04_DSP/Equalizer/FFE_DFE.m diff --git a/Classes/04_DSP/FFE_FFDCAVG.m b/Classes/04_DSP/Equalizer/FFE_FFDCAVG.m similarity index 100% rename from Classes/04_DSP/FFE_FFDCAVG.m rename to Classes/04_DSP/Equalizer/FFE_FFDCAVG.m diff --git a/Classes/04_DSP/FFE_adaptive_decision.m b/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m similarity index 100% rename from Classes/04_DSP/FFE_adaptive_decision.m rename to Classes/04_DSP/Equalizer/FFE_adaptive_decision.m diff --git a/Classes/04_DSP/VNLE.m b/Classes/04_DSP/Equalizer/VNLE.m similarity index 100% rename from Classes/04_DSP/VNLE.m rename to Classes/04_DSP/Equalizer/VNLE.m diff --git a/Datatypes/delay_mode.m b/Datatypes/delay_mode.m new file mode 100644 index 0000000..fb1736f --- /dev/null +++ b/Datatypes/delay_mode.m @@ -0,0 +1,8 @@ +classdef delay_mode < int32 + + enumeration + samples (1) + time (2) + end + +end diff --git a/Functions/calc_ber.m b/Functions/calc_ber.m index 8322d88..142ac1a 100644 --- a/Functions/calc_ber.m +++ b/Functions/calc_ber.m @@ -10,6 +10,8 @@ end options.skip_end = abs(options.skip_end); options.skip_front = abs(options.skip_front); +assert((options.skip_end+options.skip_front)