Files
imdd_silas/Classes/00_signals/Signal.m
2024-09-02 09:00:41 +02:00

754 lines
26 KiB
Matlab

classdef Signal
%SIGNAL Summary of this class goes here
% Detailed explanation goes here
properties
signal
logbook
fs
end
methods
function obj = Signal(signal,options)
%SIGNAL Construct an instance of this class
% Detailed explanation goes here
arguments
signal
options.fs = [];
end
obj.signal = signal;
obj.signal = obj.signal;
obj.fs = options.fs;
SignalType = [];
TimeStamp = [];
Length = [];
SignalPower = [];
Nase = [];
Description = [];
obj.logbook = table(SignalType,TimeStamp,Length,SignalPower,Nase,Description);
end
%% CONVERT TO INFORMATIONSIGNAL
function [i_sig, varargout] = Informationsignal(obj,options)
arguments
obj
options.fs
options.logbook
end
if isa(obj,'Electricalsignal')
%convert to information
varargout{1} = obj.fs;
i_sig = Informationsignal(obj.signal,"fs",options.fs,"logbook",options.logbook);
elseif isa(obj,'Opticalsignal')
error("Cannot convert from optical- to informationsignal. Use O/E conversion first.");
end
end
%% CONVERT TO Electricalsignal
function [e_sig, varargout] = Electricalsignal(obj,options)
arguments
obj
options.fs
options.logbook
end
obj.logbook = options.logbook;
if isa(obj,'Opticalsignal')
%convert to electrical
varargout{1} = obj.nase;
varargout{2} = obj.lambda;
e_sig = Electricalsignal(obj.signal,"fs",obj.fs,"logbook",obj.logbook);
elseif isa(obj,'Informationsignal')
try
% specify fs at varargin{1}
e_sig = Electricalsignal(obj.signal,"fs",options.fs,"logbook",options.logbook);
catch
error("Signal Conversion failed [I -> E] ");
end
end
end
%% CONVERT TO Opticalsignal
function o_sig = Opticalsignal(obj, options)
arguments
obj
options.fs
options.logbook
options.nase
options.lambda
end
fn = fieldnames(options);
for l = 1:numel(fn)
try
obj.(fn{l}) = options.(fn{l});
end
end
if isa(obj,'Electricalsignal')
%convert to optical
o_sig = Opticalsignal(obj.signal,"fs",obj.fs,"lambda",options.lambda,"logbook",obj.logbook,"nase",options.nase);
elseif isa(obj,'Informationsignal')
error("Cannot convert from information- to opticalsignal. Use E/O conversion first.");
end
end
function plot(obj, options)
% signal to plot: obj.signal
% fsamp : obj.fs (e.g. 92e9 => 92 GHz)
% length: length(obj.signal)
arguments
obj
options.fignum
options.displayname = [];
options.timeframe = 0;
end
figure(options.fignum); % If figure does not exist, create new figure
% 2) Plot into the figure handle found or created in one
t = (0:length(obj.signal)-1) / obj.fs; % time vector
if options.timeframe ~= 0
%only show a certain timeframe of signal
t = t(t<options.timeframe);
end
% 2 a) Actual plot (hold on, displayname??)
dn = options.displayname;
if isa(obj,'Opticalsignal')
sig = abs(obj.signal).^2;
else
sig = obj.signal;
end
hold on;
plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
% 2 c)
% - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)
% - ylabel amplitude
if isempty(get(gca, 'XLabel').String)
xlabel('Time (mu s)');
end
if isempty(get(gca, 'YLabel').String)
ylabel('Amplitude');
end
% Convert time axis to milliseconds for readability
xticks = get(gca, 'XTick');
set(gca, 'XTick', xticks, 'XTickLabel', xticks * 1e6);
% Add legend if not already present
if isempty(get(gca, 'Legend'))
legend;
end
hold off;
end
%% Add signals from one signal to another, the first object will sustain
function Sum = plus(X,y)
if isa(y,'Signal')
Sum = X;
Sum.signal = X.signal + y.signal;
elseif isnumeric(y)
Sum = X;
Sum.signal = X.signal + y;
end
end
function Product = times(X,y)
if isa(y,'Signal')
Product = X;
Product.signal = X.signal .* y.signal;
elseif isnumeric(y)
Product = X;
Product.signal = X.signal .* y;
end
end
%% Display length
function return_length = length(obj)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
return_length = length(obj.signal);
end
%% Write Logbook Entry
function obj = logbookentry(obj,varargin)
if nargin > 1
Description = varargin{1};
else
Description = "";
end
SignalType = [string(class(obj))];
TimeStamp = [(datetime('now','TimeZone','local','Format','HH:mm:ss'))];
Length = num2str(obj.length, ['%' sprintf('.%df', 0)]);%[obj.length];
SignalPower = [obj.power];
Nase = [0];
cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, Description};
obj.logbook = [obj.logbook;cell];
end
%% Resample Signal
function obj = resample(obj,options)
arguments
obj Signal
options.fs_in double
options.fs_out double
options.n double = 10;
options.beta double = 5;
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' ];
obj = obj.logbookentry(desc);
obj.fs = options.fs_out;
end
%%
function spectrum(obj,options)
arguments
obj
options.fignum
options.displayname = "";
end
% 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"
figure(options.fignum); % If figure does not exist, create new figure
hold on
plot(w.*1e-9,p_dbm,'DisplayName',options.displayname,'LineWidth',1);
xlabel("Frequency in GHz");
%ylabel("Power/frequency (dB/Hz)");
ylabel("Power (dBm)");
xlim([-obj.fs/2 obj.fs/2].*1e-9)
edgetick = 2^(nextpow2(obj.fs*1e-9));
% xticks([-edgetick:16:edgetick]);
xlim([-244, 244])
ylim([-120,-0]);
yticks([-200:10:10]);
legend
end
%% Power of signal
function pow = power(obj,options)
arguments
obj
options.unit power_notation = power_notation.dBm
end
pow = mean(abs(obj.signal).^2);
switch options.unit
case power_notation.dBm
if isa(obj,'Electricalsignal')
pow = pow / 50;
end
pow = 10*log10(pow)+30; %dbm
case power_notation.mW
pow = pow .* 1e3; %mW
case power_notation.W
%pow = pow % Watt
end
end
%% Peak Power of Signal
function pow_pk = power_peak(obj,options)
arguments
obj
options.unit power_notation = power_notation.dBm
end
pow_pk = max(abs(obj.signal).^2); % dBm
switch options.unit
case power_notation.dBm
pow_pk = pow2db(pow_pk)+30; %dbm
case power_notation.mW
pow_pk = pow_pk .* 1e3; %mW
case power_notation.W
pow_pk = pow_pk; % Watt
end
end
%% PAPR of signal
function papr = papr_lin(obj)
%PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power)
% divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor.
% papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2;
papr = obj.power_peak("unit",power_notation.W) / obj.power("unit",power_notation.W); %linear
end
%% PAPR of signal
function papr_db = papr_db(obj)
%PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power)
% divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor.
% papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2;
papr = obj.power_peak("unit",power_notation.W) / obj.power("unit",power_notation.W);
papr_db = 10*log10(papr);
average_power = obj.power;
peak_power = obj.power_peak;
papr_db = peak_power - average_power; %db
end
%%
function obj = normalize(obj,options)
arguments
obj Signal
options.mode normalization_mode = normalization_mode.rms
end
switch options.mode
case normalization_mode.rms
obj.signal = obj.signal/sqrt(mean(abs(obj.signal).^2,"all"));
case normalization_mode.oneone
obj.signal = obj.signal - min(obj.signal);
obj.signal = obj.signal/max(abs(obj.signal));
obj.signal = (2*obj.signal) - 1;
end
end
%%
function [obj,delay_n] = delay(obj,options)
arguments
obj Signal
options.delay_samples = 0
end
obj.signal=delayseq(obj.signal,options.delay_samples);
end
%%
function [obj,D,cuts] = tsynch(obj,options)
% time sync and cut
arguments
obj Signal
options.reference Signal
options.fs_ref = 0;
end
%normalize the signal
a = obj.normalize("mode","oneone").signal;
%resample the reference
q = obj.fs/options.fs_ref;
b = options.reference.resample("fs_in",options.fs_ref,"fs_out",obj.fs).normalize("mode","oneone").signal;
%estimate delay between signals
[co,lags] = xcorr(a,b,100);
[~,pos] = max(co);
D = lags(pos);
if D == 100
warning('Check Sync: Delay is found to be 100 which is the max. windowlength is xcorr function!')
end
% delay by lagging samples
obj = obj.delay("delay_samples",-D);
cuts = obj.length-(options.reference.length*q);
if cuts > 10
warning('Check Sync: Signal difference larger than 10.')
end
if cuts < 0
% warning('Check Sync: Reference Signal shorter than signal to sync.')
else
% then cut out the length of the reference
obj.signal = obj.signal(1:end-cuts);
end
end
function er = extinctionratio(obj,fsym,M)
histpoints = 1024; %% 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
hist_data=zeros(histpoints,histpoints_horizontal ); %% initilize eye diagram
if isa(obj,'Opticalsignal')
sig = abs(obj.signal).^2;
elseif isa(obj,'Electricalsignal')
sig = obj.signal;
else
sig = obj.signal;
end
x = (sig); %% make input signal rea)l
x = resample(x,fsym*histpoints_horizontal/2,obj.fs); %% up sample to original fsym rate
if mod(length(x),2)==1 %% if the signal lenght is not divisible by 2 (symbols displayed in the eye diagram are 2) remove last symbol
x = x(1:end-1);
end
eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols)
maxA = max(sig(100:end-100));
minA = min(sig(100:end-100));
difference= maxA-minA;
data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1;
for n=1:size(data_ind_y,1)
nn=histcounts(data_ind_y(n,:),1:histpoints+1);
hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-(
end
plot_data = 20*log10(hist_data);
plot_data(plot_data==-Inf) = 0;
maxall = 0;
for l = 1:size(plot_data,2)
[maxpk_,pos_] = max(plot_data(:,l));
if maxpk_ > maxall
maxall = maxpk_;
posxall = l;
posyall = pos_;
end
end
hist_interest = plot_data(:,posxall);
hist_interest_smoth = smooth(hist_interest,20);
[pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);
for i = 1:numel(loc)
ppeak(i) = maxA - (difference/histpoints*loc(i));
end
if isa(obj,'Opticalsignal')
er=10*log10(ppeak(1)/ppeak(end));
elseif isa(obj,'Electricalsignal')
if mean([ppeak(1),ppeak(end)]) < 1e-2
disp("No Extiction Ration for Bipolar Electrical Signal. Calculating Outer OMA instead...")
er=max(ppeak)-min(ppeak);
else
er=10*log10(ppeak(1)/ppeak(end));
end
else
er=10*log10(ppeak(1)/ppeak(end));
end
if 0
findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);
end
end
function eye(obj,fsym,M)
mode = 1;
histpoints = 1024; %% 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
hist_data=zeros(histpoints,histpoints_horizontal ); %% initilize eye diagram
if isa(obj,'Opticalsignal')
sig = abs(obj.signal).^2;
elseif isa(obj,'Electricalsignal')
sig = obj.signal;
else
sig = obj.signal;
end
x = (sig); %% make input signal rea)l
x = resample(x,fsym*histpoints_horizontal/2,obj.fs); %% up sample to original fsym rate
if mod(length(x),2)==1 %% if the signal lenght is not divisible by 2 (symbols displayed in the eye diagram are 2) remove last symbol
x = x(1:end-1);
end
eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols)
figure(922)
clf
if mode == 2
% generate "intuitive eye diagram" by drawing lines on top over
% each other; only draw 1000 lines, otherwise the plot is too
% crowded
col = cbrewer2('Set1',2);
for n=1:1000
hold on
plot(eye_mat(:,n),'LineStyle',':','LineWidth',0.1,'Color',col(2,:));
end
ylabel('Amplitude of Signal');
elseif mode == 1
% generate eye diagram using histogram
maxA = max(sig(100:end-100));
minA = min(sig(100:end-100));
% maxA = 0.0015;
% minA = 0;
difference= maxA-minA;
data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1;
for n=1:size(data_ind_y,1)
nn=histcounts(data_ind_y(n,:),1:histpoints+1);
hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-(
end
plot_data = 20*log10(hist_data);
plot_data(plot_data==-Inf) = 0;
imagesc(plot_data);
% beautify
colormap(cbrewer2("Blues",4096));
if isa(obj,'Opticalsignal')
title("Optical Eye")
ylabel("Power in mW");
y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,16));
min_ = min(abs(obj.signal(100:end-100)).^2);
max_ = abs(max(obj.signal(100:end-100)).^2);
elseif isa(obj,'Electricalsignal')
title("Electrical Eye")
ylabel("Voltage in V");
y_tickstring = string(linspace(maxA,minA,16));
min_ = min(obj.signal(100:end-100));
max_ = abs(max(obj.signal(100:end-100)));
else
title("Digital Eye")
ylabel("Digital Signal Amplitude");
y_tickstring = string(linspace(maxA,minA,16));
min_ = min(obj.signal(100:end-100));
max_ = abs(max(obj.signal(100:end-100)));
end
xlabel('Time in ps')
% add information
if 1
pwr_dbm = round(obj.power,3);
pwr_lin = obj.power("unit",power_notation.W);
papr_ = obj.papr_lin;%round(papr(obj.signal.^2),3);
yline( histpoints-(pwr_lin - minA)/difference*histpoints );
yline( histpoints-(min_ - minA)/difference*histpoints );
yline( histpoints-(max_ - minA)/difference*histpoints );
maxall = 0;
for l = 1:size(plot_data,2)
[maxpk_,pos_] = max(plot_data(:,l));
if maxpk_ > maxall
maxall = maxpk_;
posxall = l;
posyall = pos_;
end
end
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]
boxColor = [0.9 0.9 0.9]; % Light grey background color
boxEdgeColor = 'k'; % Black edge color
boxLineStyle = '--'; % Dashed line style
boxFontWeight = 'bold'; % Bold font
% Create first annotation box for Power
annotation('textbox', boxPosition, ...
'String', ['Power: ',num2str(pwr_dbm),' dBm'], ...
'BackgroundColor', boxColor, ...
'EdgeColor', boxEdgeColor, ...
'LineStyle', boxLineStyle, ...
'FontWeight', boxFontWeight, ...
'HorizontalAlignment', 'center');
% Adjust position for the second box (slightly to the right)
boxPosition = [0.37 0.86 0.2 0.05]; % Adjusted position
% Create second annotation box for PAPR
annotation('textbox', boxPosition, ...
'String', ['PAPR(lin):',num2str(papr_),''], ...
'BackgroundColor', boxColor, ...
'EdgeColor', boxEdgeColor, ...
'LineStyle', boxLineStyle, ...
'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
% 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');
yticks(linspace(0,histpoints,16));
y_tickstring = sprintfc('%.2f', y_tickstring);
yticklabels(y_tickstring);
xticks(linspace(0,histpoints_horizontal,8))
x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12);
xticklabels(x_tickstring);
end
end
% disp('h');
%
% fsig = obj.fs;
% q = fsig/fsym;
%
% if q > 10 && isinteger(q)
% sig = (obj.signal);
% else
% sig = (obj.resample("fs_in",fsig,"fs_out",fsym*30).signal);
% q = 10;
% end
%
% figure()
% clf
% cursor = 200*q;
%
% for s = 1:700
% plot(sig(cursor-q:cursor+q),'Color','black','LineWidth',0.1,'LineStyle','-');
% hold on
% cursor=cursor+q;
% s=s+1;
% end
%
% ylim([-3 3]);
end
end
end