added my matlab solutions as well as the official solutions

code/week02/Ex1_yhofmann.m

@@ -0,0 +1,46 @@
+v = randn(100,1)*sqrt(0.3) + 0.5;
+f = @(theta) -sum((y - (theta(1)*x + theta(2)*(2*x.^2 -1) + v)).^2)
+
+theta_initial = [0,0]
+options = optimset('Display', 'iter', 'TolFun', 1e-6);
+estimated_theta = fminunc(f, theta_initial, options);
+
+fprintf('Estimated theta1: %.4f\n', estimated_theta(1));
+fprintf('Estimated theta2: %.4f\n', estimated_theta(2));
+%%
+
+% Number 1
+
+for i = 1:size(x,1)
+w_k = [x(i),  2*x(i)^(2)-1];
+w(i,1) = w_k(1);
+w(i,2) = w_k(2);
+end
+
+
+%better 
+w = [x 2*x.^2-1];
+
+
+
+
+theta = inv(transpose(w)*w)*transpose(w)*(y-0.5)
+
+
+% Number 2
+var_v = 0.3;
+mu_v = 0.5;
+var_theta =0.02;
+mu_theta = [1.3; 0.9];
+theta_MAP = inv((1/var_v)*transpose(w)*w + (1/var_theta)*eye(2)) * ((1/var_v)*transpose(w)*(y-mu_v)+mu_theta/var_theta)
+
+%Number 3
+error_ML = 0;
+error_MAP =0;
+for i = 1:size(x_v,1)
+error_ML = error_ML + (theta(1)*x_v(i) + theta(2)*(2*(x_v(i))^2-1) + 0.5 - y_v(i))^2;
+error_MAP = error_MAP + (theta_MAP(1)*x_v(i) + theta_MAP(2)*(2*(x_v(i))^2-1) + 0.5 - y_v(i))^2;
+
+end
+error_ML =error_ML/size(x_v,1)
+error_MAP = error_MAP/size(x_v,1)

code/week03/Ex2_yhofmann.m

@@ -0,0 +1,110 @@
+
+clear all
+close all
+clc
+
+
+
+%% Ex 3
+close all
+
+data1 = importdata("experiment1.dat").data;
+t1 = data1(:,1),
+y1 =  data1(:,2)
+std1= data1(:,3)
+reg = [ones(length(y1),1), t1, t1.^2 ];
+theta = reg\y1
+
+
+R = diag(std1.^2);
+Z_t_direct = theta /y1; % --> Not good don't do this as probably some numerical stability issues
+% or with formula
+Z_t_formula = inv(reg'*inv(R)*reg)*reg'*inv(R);
+
+
+
+% dont do cov_theta1 = Z_t_direct*R*transpose(Z_t_direct)
+
+% Both work well
+cov_theta2 = Z_t_formula*R*Z_t_formula'
+%or direct formula
+cov_theta3 =inv(reg'*inv(R)*reg);
+
+cov1 = sqrt(cov_theta3(1,1))
+cov2 = sqrt(cov_theta3(2,2))
+cov3 = 2*sqrt(cov_theta3(3,3))
+
+figure;
+errorbar(t1,y1,std1, 'LineStyle','none','Marker','.')
+hold on
+t_cont = t1(1):0.01:t1(end);
+
+y_pred= reg*theta;
+y0 = theta(1,1)
+v0 = theta(2,1)
+a = theta(3,1)
+plot(t_cont, polyval([a v0 y0],t_cont))
+
+g_exp1 = 2*a
+
+hold off
+
+%%
+%%Ex 3 part 2
+data2 = importdata("experiment2.dat").data;
+t2 = data2(:,1),
+y2 =  data2(:,2)
+std2= data2(:,3)
+reg2 = [ones(length(y2),1), t2, t2.^2 ];
+%w = cov with formula from ex1
+R2 = diag(std2.^2);
+cov =inv(reg2'*inv(R2)*reg2);
+cov12 = sqrt(cov(1,1))
+cov22 = sqrt(cov(2,2))
+cov32 = 2*sqrt(cov(3,3))
+
+w = R2;
+
+
+%or for weighted:
+theta2 = inv(reg2'*inv(w)*reg2)*reg2'*inv(w)*y2
+%
+
+
+
+figure;
+errorbar(t2,y2,std2, 'LineStyle','none','Marker','.')
+hold on
+
+
+
+y02 = theta2(1,1)
+v02 = theta2(2,1)
+a2 = theta2(3,1)
+t_cont = t2(1):0.01:t2(end);
+plot(t_cont, polyval([a2 v02 y02],t_cont))
+
+g_exp2 = 2*a2
+
+hold off
+
+%% Ex 4
+w0 = 10;
+theta0 = pi/8;
+data3 = importdata("experiment3.dat").data;
+t3 = data3(:,1);
+y3 =  data3(:,2);
+
+figure;
+plot(t3,y3)
+
+std3= data2(:,3);
+reg3 = [ones(size(t3,1),1), t3, 0.5* t3.^2, cos(theta0 + w0*t3)];
+theta3 = reg3\y3;
+
+disp(theta3')
+
+true_val = [12.42, 44.19, -6.42, 1.53];
+disp("True values")
+disp(true_val)
+

code/week04/Ex3_yhofmann.m

@@ -0,0 +1,67 @@
+close all
+clear all
+clc
+%% Test 1
+test_Signal = randn(1,100);
+res_fft = fft(test_Signal,100);
+res_self= fourierTransformEx3(test_Signal,100);
+if(res_self ~= res_fft)
+    disp('False')
+else
+    disp('Is correct')
+end
+
+dif = sum(abs(res_self-res_fft))
+%is correct now
+
+%%
+% a)
+e = randn(1,1024);
+mean(e);
+var(e);
+
+%%
+% b)
+res = (abs(fft(e)).^2)/1024;
+
+%[pxx,f] = periodogram(e); Do not use the built in function
+%freq = linspace(0,pi,1024); Don't use linspace
+freq = 2*pi/1024*(0:1024-1);  %Use this
+idx = find(freq>=0 & freq <pi);%Necessary!!!
+
+
+figure(1)
+loglog(freq(idx),res(idx))
+xlim([0,pi])
+
+
+%%
+% c)
+T = 1;
+z = tf('z', T);
+P = (z+0.5)/(0.5+ (z+0.5)*(z-0.5)^2);
+time_steps = (0:length(e)-1);
+w = lsim(P,e,time_steps);
+
+%%
+% d)
+res_1024 = ((abs(fft(w)).^2)/1024)';
+
+%Questions: Where now only take the ones up to pi?
+%Plot goes further to the right than solutions, is solution only axis
+%limits different or 
+
+%%
+%e)
+figure(2)
+loglog(freq(idx),res_1024(idx));
+hold on
+
+[H, w] = freqresp(P, freq);
+magnitude = squeeze(abs(H));
+loglog(w(idx),magnitude(idx).^2);
+hold on
+
+dif = res_1024' - magnitude;
+loglog(freq(idx),abs(dif(idx)));
+legend('Periodogram','Plant', 'Difference')

code/week04/fourierSeriesEx3.m

@@ -0,0 +1,24 @@
+% function f =fourierSeriesEx3(u, N)
+% 
+% f = zeros(N,1);
+% for n = 1:N
+%     w_n = 2*pi *(n-1)/N;
+%     temp = 0;
+%     for k = 1:n
+%         f(n) = temp + u(k)*exp(-1j*w_n*(k-1));
+%         temp = f(n);
+%     end
+% 
+% end
+
+% end
+
+function [X] = fourierSeriesEx3(x,N)            
+    X = zeros(1, N);            % Initialize the DFT result
+    for k = 1:N
+        X(k) = 0;
+        for n = 1:N
+            X(k) = X(k) + x(n) * exp(-1j*2*pi*(k-1)*(n-1)/N);
+        end
+    end
+end

code/week04/fourierTransformEx3.m

@@ -0,0 +1,9 @@
+function [X] = fourierTransformEx3(x,N)            
+    X = zeros(1, N);            % Initialize the DFT result
+    for k = 1:N
+        X(k) = 0;
+        for n = 1:N
+            X(k) = X(k) + x(n) * exp(-1j*2*pi*(k-1)*(n-1)/N);
+        end
+    end
+end

code/week05/Ex4_yhofmann.m

@@ -0,0 +1,123 @@
+close all
+clear all
+clc
+%%
+%a)
+N = 1024;
+e = randn(N,1)*sqrt(0.01); % pay attention is square root for the standard deviation!!!!!!!!
+u = 2 * randi([0, 1], 1, N) - 1;
+
+
+T = 1;
+z = tf('z', T);
+H = 0.5*(z-0.9)/(z-0.25)
+G = (0.1*z)/(z^4-2.2*z^3+2.42*z^2-1.87*z+0.7225)
+time_steps = T*(0:N-1);
+resu = lsim(G,u,time_steps);
+rese = lsim(H,e,time_steps);
+y_experimet = resu + rese;
+omega = (2*pi/(T*N)*(0:N-1))';
+idx = find(omega>=0 & omega< pi); 
+
+% Uncomment to see plots of signals
+%{
+figure(1)
+plot(u)
+hold on
+plot(e)
+hold on
+plot(y_experimet)
+legend('u', 'e', 'Y')
+
+
+figure(2)
+plot(rese)
+hold on
+plot(resu)
+hold on
+plot(y_experimet)
+hold on
+legend('H(e)', 'G(u)', 'Y')
+%}
+
+%%
+%b)
+U = fft(u);
+Y_exp = fft(y_experimet);
+
+% if one of them has many which are 0 is probably due to periodic input -->
+% only calculate at these
+
+ % figure(1)
+ % loglog(omega(idx), abs(U(idx)));
+% title('U')
+% figure(2)
+% loglog(omega(idx), abs(Y_exp(idx)));
+% title('Y_exp')
+
+G_est = Y_exp ./transpose(U); %Watch out as U is complex don't do U' but use transpose(U)!!!!!!!!!
+G_fresp_true = squeeze(freqresp(G,omega));
+
+
+figure(3)
+loglog(omega(idx),abs(G_est(idx)))
+hold on
+loglog(omega(idx),abs(G_fresp_true(idx)))
+hold on
+Err = G_est - G_fresp_true;
+loglog(omega(idx),abs(Err(idx)))
+legend('Estimate', 'True', 'Error');
+title('Magnitude')
+
+
+
+
+%%
+%c)
+N_c = N/4;
+omega_c = (2*pi/(T*N_c)*(0:N_c-1))';
+idx_c = find(omega_c>=0 & omega_c<= pi);
+
+% Does the same but not really good to do it like that
+Y_experimet1 = fft(y_experimet(1:end/4)); 
+Y_experimet2 = fft(y_experimet(end/4+1:end/2));
+Y_experimet3 = fft(y_experimet(end/2+1:3*end/4));
+Y_experimet4 = fft(y_experimet(3*end/4+1:end));
+U1 = fft(u(1:end/4));
+U2 = fft(u(end/4+1:end/2));
+U3 = fft(u(end/2+1:3*end/4));
+U4 = fft(u(3*end/4+1:end));
+G1_est = Y_experimet1./transpose(U1);
+G2_est = Y_experimet2./transpose(U2);
+G3_est =Y_experimet3./transpose(U3);
+G4_est =Y_experimet4./transpose(U4);
+
+%Better like that
+u_res = reshape(u,N_c,4);
+y_res = reshape(y_experimet,N_c,4);
+
+U_res = fft(u_res,[],1);
+Y_res = fft(y_res,[],1);
+
+G_est_res = Y_res./U_res;
+
+
+G_est_average = zeros(length(G1_est),1);
+for k =1:length(G1_est)
+    G_est_average(k) = mean([G1_est(k);G2_est(k);G3_est(k); G4_est(k)]);
+end
+G_sum = G1_est+G2_est+G3_est+G4_est;
+G_est_average2 = (G_sum)/4; % both ways work now
+
+
+G_est_average2 = mean(G_est_res,2);
+
+
+figure(5)
+loglog(omega_c(idx_c),abs(G_est_average2(idx_c)))
+hold on
+loglog(omega(idx),abs(G_est(idx)))
+hold on
+loglog(omega(idx),abs(G_fresp_true(idx)))
+legend('Averaged estimate', 'Estimate', 'True System');
+title('Magnitude from averaged')