[u,y,koef,lamda]=systemdesign(N,optu,optrn,por)

function [u,y,koef,lamda]=systemdesign(N,optu,optrn,por)

%y(k)=-Sum[a(i)*yi(k)]+ Sum[b(j)*uj(k)]+rn;
na=8;
nb=8;
if nargin<1, 
N=1e3;
end;
if nargin<4,
por=5;
end; 
if nargin<2,
    optu=1;
end
if optu==1
             u=linspace(-N/2,N/2)';
%            u=linspace(-pi,pi)'*cos(-pi:pi);
%            u = firls(N-1,[0 .4 .5 1],[1 1 0 0],[1 10]); 
    else if  optu==2
             u = random('Poisson',1:N,1,N);
        else if optu==3
             u = random('Normal',0,1,1,N);
             end
        end
end
 
if nargin<3,
    optrn=1;
end
if optrn==1
    rn = firls(80,[0 .4 .5 1],[1 1 0 0],[1 10]);
    else if  optrn==2
         rn = random('Poisson',1:N,1,N);
        else if optrn==3             
         rn = random('Normal',0,1,1,N);
             end
        end
end

por=por/100;
rn=rn*(1+randn(1,1)*por);

s=tf('s');
syms yt ys f0 A ex
y=zeros(8,8);  f0=zeros(8,1); A=zeros(8,1);ex=zeros(8,1);
for k=1:8,                                                              
H(k) = k/(k^3*s^2+k^2*s+k);
freq=1/N;
w=2*pi*freq;
f0(k,1)=atan((k^2-k^3*w)/k);
ex(k,1) = exp(j*(2*pi*freq*k+f0(k,1)));
A(k,1)=sqrt(k^2/[(k^2-k^3*w^2)^2+k^2]);
y=y+A*ex';
end
figure
plot(y)
xlabel('Abtastzeiten'),ylabel('Ausgangsignal [m]')
title('Ausgangsignalverlauf');

for k=1:8,
z=zeros(1,4);
a=k^3; b=k^2; c=k;
discri = b*b-4.*a*c;
if (discri >=0)
z(1)=(-b-sqrt(discri))/2./a;
z(2)=(-b+sqrt(discri))/2./a;
else
z(3)=[-b+sqrt(-1.*discri)]/2*a; 
z(4)=[-b-sqrt(-1.*discri)]/2*a;
end
end
figure
plot(z)
xlabel('Abtastzeiten'),ylabel('Eigenform [m]')
title('Polenstellenplot - Eigenformverlauf');

sig_laen  = length(u);
mittel_qdr_sig=[1,sig_laen];
rn_laen = length(rn);
mittel_qdr_sig=[1,rn_laen];
snrwert=[1,rn_laen];

for l=1:rn_laen;
mittel_qdr_sig(l)  = mean(sum(u(l)^2,2)/sig_laen); 
mittel_qdr_rn(l) = mean(sum(rn(l)^2,2)/rn_laen);
snrwert(l) = mittel_qdr_sig(l)/mittel_qdr_rn(l);
end

figure,
subplot(211), plot(mittel_qdr_sig,'b-'),grid 
xlabel('Abtastzeiten'),ylabel('Eingangsignalmittelwert [dB]')
title('Eingangsignalmittelwertverlauf');
subplot(212), plot(mittel_qdr_rn,'r-'),grid 
xlabel('Abtastzeiten'),ylabel('Rauchsignalmittelwert')
title('Rauchsignalmittelwertverlauf');
 
figure,
bode(H), grid
xlabel('Kreisfrequenz'),ylabel('Amplitude') 
title('Phasengang: Phi(jw) & Amplitudengang: G(jw)');

%----------Vergissenheitsfaktor----------
lamda= random('Normal',-1,1,1,N);
for k=1:length(lamda),
    lamda(k) = lamda(1)*lamda(k)+lamda(length(lamda))* (1-lamda(1));
end

figure
subplot(211), plot(log(abs(lamda)),'b-'),grid 
xlabel('Abtastzeiten'),ylabel('lamda [dB]')
title('lamda-verlauf');
subplot(212), plot(lamda,'r-'),grid 
xlabel('Abtastzeiten'),ylabel('lamda')
title('lamda-verlauf');

save datensatz.mat y u H rn snrwert lamda na nb

end

function [eigw,zeta,modaldaemp]=modalparameter(tabtast)
if nargin<1,   tabtast=1e-3; end
for k=1:8,
z=zeros(1,4);
a=k^2; b=k; c=1;
discri = b*b-4.*a*c;
if (discri >=0)
z(1)=(-b-sqrt(discri))/2./a;
z(2)=(-b+sqrt(discri))/2./a;
else
z(3)=[-b+sqrt(-1.*discri)]/2*a; 
z(4)=[-b-sqrt(-1.*discri)]/2*a;
end
end
for i=1:length(z)
eigw(i)=sqrt(log(z(i))*log(conj(z(i))))/tabtast;%eigenvalue
zeta(i)=-(log(z(i))+log(conj(z(i))))/(2*eigw(i)*tabtast);%madal daemping value
modaldaemp(i)=2*eigw(i)*zeta(i);%modal daemping
end

save modalparameter.mat eigw zeta modaldaemp
end