WIP büro

This commit is contained in:
Silas Oettinghaus
2023-05-26 11:06:39 +02:00
parent 8b2129ac03
commit b640c013cb
20 changed files with 192 additions and 52 deletions

162
Classes/02_etc/Amplifier.m Normal file
View File

@@ -0,0 +1,162 @@
classdef Amplifier
%UNTITLED Summary of this class goes here
% Detailed explanation goes here
properties
amp_mode
gain_mode
nase_mode
amplification_db
noifig
fsimu
end
methods
function obj = Amplifier(options)
%UNTITLED Construct an instance of this class
% Detailed explanation goes here
arguments
options.amp_mode amp_mode = amp_mode.ideal_no_noise
options.gain_mode gain_mode = gain_mode.output_power
options.nase_mode nase_mode = nase_mode.pass_ase
options.amplification_db double = 0
options.noifig double = 0;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
end
function signalclass_out = process(obj,signalclass_in)
% actual processing of the signal
signalclass_in = obj.process_(signalclass_in);
% append to logbook
lbdesc = ['Amp '];
signalclass_in = signalclass_in.logbookentry(lbdesc);
% write to output
signalclass_out = signalclass_in;
end
function [X_out] = process_(obj,X_in)
%calc gain for output power mode, amp mode and saturated mode
a_lin = obj.calculateGain(X_in.signal);
%apply amplification
X_in.signal = a_lin * X_in.signal;
% now, optical ase noise:
if class(X_in) == "Opticalsignal"
if obj.amp_mode == amp_mode.ideal_no_noise
%don't add/ remove noise, but scale it accordingly
X_in.nase = obj.onlyAmplifyAse(X_in.nase, a_lin);
elseif obj.amp_mode == amp_mode.edfa_increase_nase
%calculate ASE-noise
X_in.nase = obj.increaseAse(X_in.nase, a_lin, obj.noifig , X_in.lambda );
elseif obj.amp_mode == amp_mode.edfa_replace_nase
%calculate ASE-noise
X_in.nase = obj.replaceAse(X_in.nase, a_lin, obj.noifig , X_in.lambda );
end
if obj.nase_mode == nase_mode.generate_ase
nase = X_in.nase;
fs = X_in.fs;
dimension = size(X_in.signal);
nase_numeric = obj.generateAseNoise(nase, fs, dimension);
[X_in.signal, osnr] = obj.applyAseToSignal(X_in.signal, nase_numeric);
X_in.nase = 0;
elseif obj.nase_mode == nase_mode.pass_ase
X_in.nase = X_in.nase;
end
end
X_out = X_in;
end
function a_lin = calculateGain(obj,xin)
if obj.gain_mode == gain_mode.output_power
%get linear gain for output power mode
pow_in = mean(abs(xin.^2)) ; % lin input power
pow_out = 10^(obj.amplification_db/10 - 3) ; % dBm to lin
a_lin = sqrt(pow_out/pow_in) ;
elseif obj.gain_mode == gain_mode.gain
%get linear gain for classic gain mode
a_lin=10^(obj.amplification_db/20);
end
end
function nase_amped = onlyAmplifyAse(~,nase_old, a_lin)
nase_amped = a_lin^2 * nase_old;
end
function nase_new = calculateAseFromThisAmp(~, a_lin, noisefig, lambda)
h = Constant.Planck;
c = Constant.LightSpeed;
nase_new = 0.5* 10^(noisefig/10) * h*c/lambda * (a_lin^2-1) ;
end
function ase_increased = increaseAse(obj,nase_old, a_lin, noisefig, lambda)
nase_fromThisAmp = obj.calculateAseFromThisAmp(a_lin, noisefig, lambda);
nase_old = a_lin^2 * nase_old;
ase_increased = nase_old + nase_fromThisAmp;
end
function nase_new = replaceAse(~,a_lin, noisefig, lambda_nm)
nase_fromThisAmp = obj.calculateAseFromThisAmp(a_lin, noisefig, lambda_nm);
nase_new = nase_fromThisAmp;
end
function nase_numeric = generateAseNoise(~, nase, fs, dimension)
nase_numeric = (randn(dimension) + 1i*randn(dimension))*sqrt(nase/2*fs) ;
end
function [noisy_sig, osnr] = applyAseToSignal(~, opticalsignal, nase_numeric)
noisy_sig = opticalsignal + nase_numeric ;
osnr = pow2db(mean(abs(opticalsignal.^2))/mean(abs(nase_numeric.^2)));
disp(osnr);
end
end
end

204
Classes/02_etc/Filter.m Normal file
View File

@@ -0,0 +1,204 @@
classdef Filter
%FILTER Summary of this class goes here
% Detailed explanation goes here
properties
H
filterType
f_cutoff
lowpass
signal_length
filtdegree
passband_ripple
stopband_ripple
fsamp
w
end
methods
function obj = Filter(options)
%FILTER Construct an instance of this class
% Detailed explanation goes here
arguments
options.filterType filtertypes = filtertypes.bessel_inp ;
options.f_cutoff = 0;
options.fsamp = 0;
options.filtdegree = 3;
options.passband_ripple = 0.5;
options.stopband_ripple = 0.5;
options.lowpass = 1;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
end
function signalclass_out = process(obj,signalclass_in)
% actual processing of the signal
signalclass_in.signal = obj.process_(signalclass_in.signal);
% append to logbook
filterdesc = [num2str(obj.filtdegree),'. order ',char(obj.filterType),' filter with f_cutoff at ', num2str(obj.f_cutoff*1e-9), ' GHz.'];
signalclass_in = signalclass_in.logbookentry(filterdesc);
% write to output
signalclass_out = signalclass_in;
end
function yout = process_(obj,xin)
obj.signal_length = length(xin);
[obj.H,obj.w] = obj.buildFilter(obj.filterType);
yout = obj.applyFilter(xin);
end
function y_filtered = applyFilter(obj,xin)
y_filtered = ifft(obj.H.*fft(xin));
end
function [H,w] = buildFilter(obj,filterType)
w = [];
rp = obj.passband_ripple; %passband ripple
rs = obj.stopband_ripple; %stopband ripple
switch filterType
case 1
% Bessel filter, impulse invariant transformed
[B, A] = besself(obj.filtdegree, 2*pi*obj.f_cutoff);
[B ,A] = impinvar(B,A,obj.fsamp);
case 2
% Bessel filter, impulse bilinear transformed
[Z, P, K] = besself(obj.filtdegree, 2*pi*obj.f_cutoff);
[Z ,P, K] = bilinear(Z,P,K,obj.fsamp);
[B ,A] = zp2tf(Z ,P ,K);
case 3
% Butterworth filter
if obj.lowpass == 1 %lowpass
[B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fsamp/2),'low');
else % highpass
[B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fsamp/2),'high');
end
case 4
% Chebyshev 1 filter
[B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fsamp/2));
case 5
% Chebyshev 2 filter
[B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fsamp/2));
case 6
% Elliptic filter
[B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fsamp/2));
case 7
% Hamming filter
g=(obj.filtdegree-1)/2;
wc=obj.f_cutoff/(obj.fsamp/2);
B = wc*sinc(wc*(-g:g)).*hamming(obj.filtdegree)';
A=1;
case 8
% Raised Cosine filter
B = firrcos(obj.filtdegree,obj.f_cutoff,para.df,obj.fsamp);
A=1;
case 9
% Sinc filter
g=(obj.filtdegree-1)/2;
wc=obj.f_cutoff/(obj.fsamp/2);
B = wc*sinc(wc*(-g:g));
A=1;
case 10
% Gaussian Filter
%check if order ist multiple of 1/2
blocklen=obj.signal_length;
faxis=linspace(-obj.fsamp/2,obj.fsamp/2,blocklen+1)';%generates arow vector faxis of blocklen+1 points linearly spaced between and including -para.fs/2 and para.fs/2
faxis=ifftshift(faxis(1:end-1));
H=exp(-((faxis)/(obj.f_cutoff*2)).^(2*obj.filtdegree)*log(2)*2^(2*obj.filtdegree-1));
% figure()
% hold on
% xline(obj.f_cutoff*1e-9,'LineWidth',3,LineStyle='--');
% xline(-obj.f_cutoff*1e-9,'LineWidth',3,LineStyle='--');
% plot(faxis*1e-9,20*log10(abs(H)),'LineWidth',3);
% ax = gca;
% ylim([-6 0])
% grid on
% xlabel('Freq in GHz')
% ylabel('Magnitude (dB)')
end
% Build Filter from coefficients
if filterType ~= 10
[H,w] = freqz(B, A, obj.signal_length,'whole');
% figure()
% hold on
% plot(w/(2*pi)*obj.fsamp*1e-9,20*log10(abs(H)),'LineWidth',3);
% ax = gca;
% ylim([-6 0])
% grid on
% xlabel('Freq in GHz')
% ylabel('Magnitude (dB)')
% freqz(B, A)
%hfvt = fvtool(B,A);
end
end
function show(obj)
obj.signal_length = 1024;
[H,w] = obj.buildFilter(obj.filterType);
plot(w/pi,20*log10(abs(H)));
ax = gca;
ax.XTick = 0:.5:2;
grid on
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')
%freqz(H);
end
end
end